NASA CONTRACTOR 
REPORT 


'13-JL5Q2. 

ASA CR-2195 







EFFECT OF SURFACE ASPERITY A ^ 

ON ELASTOHYDRODYNAMIC LUBRICATION ^ 


by Kwan Lee and H. S. Cheng 


Prepared by 

NORTHWESTERN UNIVERSITY 

Evanston, III. 

for Lewis Research Center 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • FEBRUARY 1973 




1. Report No. 2. Government Accession No. 

NASA CR-2195 

3. Recipient's Catalog No. 

4. Title and Subtitle 

EFFECT OF SURFACE ASPERITY ON 
ELASTOHYDRODYNAMIC LUBRICATION 

5. Report Date 

February 1973 

6. Performing Organization Code 

7. Author(s) 

Kwan Lee and H. S. Cheng 

8. Performing Organization Report No. 

None 

10. Work Unit No. 

9. Performing Organization Name and Address 
Northwestern University 
Evanston, Illinois 

11. Contract or Grant No. 
NGL- 14-007 -084 

13. Type of Report and Period Covered 
Contractor Report 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

15. Supplementary Notes 


Project Manager, Erwin V. Zaretsky, Fluid System Components Division, NASA Lewis 
Research Center, Cleveland, Ohio 

16. Abstract 

The important aspects of elastohydrodynamic lubrication, with a single, one -dimensional 
asperity, have been found by solving numerically the coupled transient Reynolds equation and 
the elasticity equation. Even though the assumption of a single asperity is highly ideal, this 
study sheds some light on the effect of surface roughness on elastohydrodynamic lubrication. 
The results show that the film pressure tends to increase more than the steady state pressure, 
and in particular, the increase in pressure reaches a maximum as the asperity approaches the 
inlet of the contact region. The asperity height and the pressure increase above the steady 
state pressure are closely related to each other; the higher the asperity height, the larger 
the pressure increase. In the pure rolling case, it has been found that a local pressure peak 
is not developed. However, in the cases of sliding and rolling, a small, local pressure peak 
is developed on the pressure profile when the asperity moves into the contact region. In 
general, the overall film thickness profile increases with increasing asperity height, but is 
not significantly affected by the asperity width. Moreover, the slope of the overall film thick- 
ness profile for the transient cases is much greater than the steady state profile, which is 
approximately constant across the contact width. The increase in the center film thickness 
also depends upon the width and height of the asperity. Even for the case of an asperity height 
of 2H g , the center film thickness increases more than 100 percent compared to the steady 
state center film thickness. 

17. Key Words (Suggested by Author(s)) 

Elastohydrodynamics 
Lubrication 
Asperity interaction 
Surfaces 


19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified 

Unclassified 

75 

$3.00 


18. Distribution Statement 

Unclassified - unlimited 


For sale by the National Technical Information Service, Springfield, Virginia 22151 






















TABLE OF CONTENTS 


Page 

CHAPTER 1 - INTRODUCTION 

1.1 Introduction 1 

CHAPTER 2 - MATHEMATICAL FORMULATION 

2.1 Introduction 3 

2.2 Geometry of Asperity 4 

2.3 Governing Equations 4 

2.3.1 Film Thickness 4 

2.3.2 Hydrodynamic Equation. 6 

2.4 Viscosity and Density Variations 10 

2.5 Formulation of Elas tohydrodynamic 

Problem 10 

2.5.1 Coupled Time-Dependent 

Elas tohydrodynamic Equations 10 

2.5.2 Normalization 12 

2.6 Method of Solution 14 

2.6.1 Outline of Approach. 14 

2.6.2 Integration of Pressure in the 

Inlet Pressure 15 

2.6.3 Calculation of Deformation ...... 16 

2.6.4 El as tohydrodynamic Equation in 

the Middle Region 18 

2.6.5 Calculation of Center Film 

Thickness. ...... . 21 

2.6.6 Outline of Numerical Procedure .... 22 

CHAPTER 3 - DISCUSSION OF RESULTS 

3.1 Introduction 25 

3.2 Pressure Profile 25 

3.3 Film Thickness 27 

3.4 Load. ...... 29 

CHAPTER 4 - SUMMARY OF RESULTS . 30 

APPENDIX A - Symbols List 45 

APPENDIX B - Calculation of Matrix Elements in Eq. (51) .... 50 

APPENDIX C - Computer program Flow Diagram and Fortran 

Listings 55 

REFERENCES 74 


iii 



CHAPTER I 


INTRODUCTION 


1.1 INTRODUCTION 

Since Reynold developed the hydrodynamic lubrications theory, bear- 
ing performance between two parallel or two conformable surfaces can be 
readily determined by solving the Reynolds equation for the pressure 
distribution within the lubricant film. However, for highly loaded, 
counterf ormed contacts such as ball, roller, or gear tooth contacts, 
it was found that hydrodynamic theory alone cannot explain the lubrica- 
tion phenomenon. The film thickness obtained from the hydrodynamic 
equation is so small that direct metallic contacts between asperities 
must take place and cause surface distress. However, in practice, 
these concentrated contacts do operate quite satisfactorily without any 
signs of distress in lubrication. This indicates strongly that some 
degrees of hydrodynamic lubrication must exist between these counter- 
formal contacts. The analytical proof of the surface separation be- 
tween these contacts was first given by Grubin [ 1 3 by including the 
surface deformation and the effect of a variable viscosity. His work 
opened a new branch of lubrication known as elastohydrodynamic lubri- 
cation. 

During the past two decades, rapid progress has been made in elasto- 
hydrodynamic lubrication. There now exist several numerical solutions to 
determine the pressure and film thickness profiles [2 to 5], Also 
available are simplified formulas to compute the minimum film thickness 
in terms of load, speed, and lubricant parameters. However, in all 
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these EHD theories, the contacting surfaces are assumed to be perfectly 
smooth. In reality, surfaces are never perfectly smooth, and for con- 
ditions where the average height of the asperities approaches or ex- 
ceeds the mean film thickness, the smooth-film elastohydrodynamic 
theories are no longer valid. The asperities can have a significant 
influence upon the pressure and film thickness profiles. To determine 
the interaction between the lubricant and asperities, it is necessary 
to solve the transient elastohydrodynamic equations as the asperities 
moving through the contact region. 

The present study concentrates on the interaction of a single as- 
perity with lubricant as it enters the inlet region of an elastohydro- 
dynamic contact. By solving the coupled elasticity and hydrodynamic 
equations at successive time intervals during the entrance of a single 
asperity, the modification of the pressure distribution and the film 
thickness level around the asperity can be determined. 
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CHAPTER II 


MATHEMATICAL FORMULATION 


2.1 Introduction 

Due to the presence of the asperity on the moving surface, both 
the pressure and film deformation become time-dependent as the asperity 
enters the contact region. To determine the change in pressure and de- 
formation caused by the asperity, it is required to solve the coupled 
elastohydrodynamic equation at successive time intervals taking into 
account the effect of the squeeze film term. For each time interval, 
the effort required for determining the pressure and deformation pro- 
file is equivalent to a single, conventional, elastohydrodynamic solu- 
tion. 

For heavily loaded contacts, the conventional EHD theories show 
that the pressure and deformation profiles in the central region of the 
contact are almost identical to the dry-contact, Hertzian profiles. 
Deviations from Hertzian distributions only occur at the inlet and exit 
regions. This fact enables one to investigate the effect of asperity 
in the inlet and exit regions separately. It is assumed that the dis- 
turbance caused by the asperity at the inlet region is not felt in the 
exit region and vice versa. 

The present study is primarily concerned with the effect of the 
asperity in the inlet region. In solving the time-dependent EHD equa- 
tions, the pressure distribution in the exit-half of the contact region 
is assumed to be the Hertzian elliptical distribution. 
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2.2 Geometry of Asperity 


The surface geometry adopted in the present study is the per- 
fectly smooth contact surface of the cylinder attached with a single, 
one-dimensional asperity of parabolic shape as shown in Fig, 1 . In 
the present study, the height and the width of asperity are changed in 
such a way that each effect can be investigated separately„ The max- 
imum values of height and width of the asperity are twice the center 
film thickness and one-half of the Hertzian width of the same center 
pressure. The average of the rolling speeds of the two cylinders is 
kept constant which in turn fixes the center film thickness of the 
steady state solution with a fixed center pressure. 

2.3 Governing Equations 

2.3.1 Film Thickness 

The goemetrical configuration of two cylinders can be shown 
to be equivalent to a cylinder and a flat surface as shown on 
Fig. 1(c). The radius of the equivalent cylinder is 


R = 


R 1 + R 2 
R 1 R 2 


( 1 ) 


Since the contact region is much smaller compared to the radius 
of the cylinder, the geometrical film thickness without the height 
of the asperity is. 


h 

g 


= h' 


o 



( 2 ) 
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The height of the asperity is approximated by the parabolic equa- 
tion as 

f = f ^cos 9 = [(x 3 - x) 2 - c 4 2 J cos 0 (3) 

Because of the narrowness of the contact width, cos 9 = 1.0. 

It follows that the asperity height can be approximated by. 



Thus, the geometrical film thickness with the asperity height is, 



The deformation of the contact surfaces, as derived in Ref. [6], 


<= I 

d(x,t) = - ^ ^ P(§,t)£ n ' X * d£ 


(6) 


where 

A . I (— i! + — J l) 

E 2 \ E 1 E 2 ' 

and E^, and v^, are the Young's modulus and the Poisson's 
ratio of the cylinder 1 and 2, respectively. 

The film thickness including the deformation becomes 


h(x,t) = h Q (t) + 2r nE 




-CO 


nr 


di 


i r. 2 2 i 

+ 2? L (x 3" x) ‘ °4 J 


= h 1 (x,t) + f . 


(7) 
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2.3.2 Hydrodynamic Equation 

In the present study, the transient, one dimensional 
Reynolds equation is taken as the governing equation for the 
pressure distribution. 

3 ( Ph 3 3p\ u l + u 2 3(ph) . 3(ph) 

3^ Vl2 ^ Sx/‘ 2 3x + 3t (8 

As mentioned in Section 2.1, the present study is mainly 
concerned with the first half of the contact region, and the 
pressure distribution in the second half of the contact region 
is assumed to be Hertzian. Thus, the appropriate boundary con- 
dition is that the center pressure gradient is zero, but the 
center film thickness is allowed to change as dictated by the 
pressure distribution in the first half of the contact region. 

The boundary conditions for Eq. (8) are 


P = 0 


ap 

3x 


0 


at x = -“ 
at x = 0 


(9) 


(a) The Case For Pure Rolling 

When the two cylinders are in rolling motion without sliding, 
the time derivative and the spatial derivative of "f" cancel out 
each other as shown below: 


3 (ph) _ | 

3t 3t L 


P(h 

o 


3(ph Q ) 

3t 


+ f 


2R 


+ d + 


f )J 


S_p , u 

3t r 



x) 


( 10 ) 
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Since the cylinder surfaces are approximately parallel in the 
contact region, the sum of the cylinder curvature and deformation 
terms is very small. Furthermore, their time derivative is neg- 


ligible. 


U l + U 2 d(ph) 


dx 


= u 


dx 


P(hj+ f)j 


d(ph ) 

1 , , dp u s 

- u —ST - + u£ S - 7 p(x 3- x) 


In the above manipulations, the following are used: 
^ x 3 

dT = U ! and u i = u 2 = u 


(ID 


the sum of Eq. (10) and (11) is 


( 12 ) 


Substituting Eq. (12) into Eq. (8), one obtains 


d_ ( Ph 3 dP^ _ (ph l ) 
dx \L2p, dx > U dv 


a<Ph,) a(Ph o ) a p a P 

“St + £< " a! + a^ 


(13) 


where 


1 r 2 2 ”i 

*(*>* 3 ) 0 2? L (X 3‘ x) - c 4 j 


|x 3 - x| ^ c 4 


f(x,x 3 ) = 0 


l X 3' x l > c 4 


Eq. (13) is integrated from x = -x to x = 0 using the second 
boundary condition of (9). 


ap i2n ; / , ,1 \ ]° r 8<ph o > 

3x 3 L u \ ph l' ph l ) ~ .1 L at 

ph 'x=o -x 


+ £ <"|j + i>J dx } (W > 
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A new variable Q is introduced to calculate the inlet pres- 
sure and the center film thickness. The advantage of using Q is 
to eliminate strong dependence of pressure on viscosity in the 
governing equation. The explicit viscosity term does not appear 
in the governing equation if the pressure gradient is replaced 
by the derivative of Q. 

Defining Q as 


Q = l 


1 

M- 


= and 
^s 

The spatial derivative of Q is 


where p = - t — and p is the ambient viscosity. 
Lb S 

S 


(15) 


SQ _ 1 S(ln p) SP 

Sx — SP Sx 
M- 

(16) 

_ H dP 
” — Sx 
V> 


The pressure derivative in Eq. (14) is replaced by 





12p a 
s 




Ph, 


I ) " J 

x=o -X 


r° r s (ph G ) 


St 


+ f(u 


SP 

Sx 



Integrating Eq. (17) again from x = -°° to x gives 
Q - 12m. or f {(-ij) [»(ph - Ph I ) - f _■ 

_ m I 


Ph' 


x=o 


S(P V , rf Sp 

St + f(u s§ + 



(18) 


The above equation will be used in the calculation of the inlet 
pressure and the center film thickness. 
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The instantaneous load per unit width of the cylinder is the 
sum of the integrals of the first half pressure profile and the 


Hertzian pressure in the second half of the contact region. 

2 

r TTRP 
P(x,t)dx + 0 


W(t) 


(19) 


_CO 


(b) The Case For Rolling And Sliding 

When the two cylinders undergo a rolling and sliding motion, 
the time derivative and the spatial derivative of "f" do not cancel 
out each other as was the case for pure rolling. Thus, the as- 
perity effect on the fluid film becomes much stronger than that 
of the pure rolling. 

The expanded forms of the rolling and squeezing terms are: 

+£ 3£ + !i ,, 

^ ~ + f ** + - (> 


at 

u,+ u 


at 


r x) 


( 20 ) 


r u 2 a(ph) _ U 1 + u 2 r^py ai£fi 

2 dx 2 L dx 


9x J 




( 21 ) 


In Eq. (20) the sum of cylinder curvature and deformation 
terms is neglected as explained before. Substituting Eq. (20) 
and (21) into Eq. (8), one obtains 


3_ (W ap'i , u i* u 2 

ax 'i2(i ax/ 2 l ax 


+ f 


api 

ax J 


+ ^ + £ |e 


+ (u x - u 2 )(x 3 - x) 


( 22 ) 


where 
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f (x 3 ,x) = ^ [( x 3 - x ) 2 - C 4 J 


' x 3 ~ X ' ^ C 4 


f(x 3 ,x) = 0 


( V u 2 >(x 3 ' x> 2 ? 


= 0 


l x 3 " x l > C 4 


Eq. (22) is integrated from x = -x to x = 0 using the second 
boundary condition of (10). Thus we obtain 
+ u„ / \ r*° 1-11,+ u 


ap i2n TV “2 ( , \ r°rV u 2 3 p 

si ‘ T5 in — l p V ph i >- J l~ 2 (£ ai> 

ph x=o -x 


S(ph o ) 

~ 


3h 

+ f + (u, — u 0 )(x - x) 


fc]-} 


(23) 


2.4 Viscosity and Density Variations 

The property of lubricant is assumed to depend upon pressure only; 
the viscosity is a exponential function of pressure with a suitable 
pressure-viscosity coefficient and the density function used in the 
present investigation is the one developed by Dowson and Whitaker [5]. 

cyp 

H = (24) 

p - P S ( l + TT^p) < 25 > 

where p and p are the ambient viscosity and density, respectively, 
s s 


2.5 Formulation of Elastohydrodynamic Problem 

2.5.1 Coupled Time-Dependent Elastohydrodynamic Equations 

No previous work has ever undertaken a time-dependent EHD 
problem of pure rolling or rolling and sliding. The complexity 
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of the present problem is multiplied by bringing in the asperity 
action on the fluid film in the contact region. The successful 
attempt for obtaining the full solution requires not only solving 
the hydrodynamic equation and the elasticity equation, but also 
the center film thickness calculation at each instantaneous loca- 
tion of the asperity. 

The major equations to be solved are: 

(a) The case for pure rolling 


S - l u ( ph r *1 

Ph 


s *&]*■} cu, 


x=o -X 


x 2 4 r F _ x I 

h(x,t) - hjt) +|^ - ^ J P(5,t)An dC + f(x,x 3 ) (7) 

_ ao ' 


where 


f (x,x 3 ) = ~ [(x 3 - X) 2 - c J 
f(x,x 3 ) = 0 



X c 
X !> c 


4 

4 


(b) The case for rolling and sliding 



rV "2 


< £ 


5(Ph Q ) 

~Tt~ 


( 22 ) 


Z pCO 

h (x,t) = h Q (t) + — - ^ j P(?,t)£n + f (x,x 3 ) (7) 

where 
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f (x,x 3 ) 



f(x,x 3 ) = 0 

<V U 2 )(X 3‘ X) fe = ° 



> C, 


2.5.2 Normalization 

The following non-dimensional variables are introduced: 



where a is the Hertzian half-width and the subscript "o" indicates 
the variables at the film center. 


The normalized governing equation are written as: 
(a) The case for pure rolling 
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48Ua 
Q = — 


rOC 


H H 


m*¥)-=±r £&+■*<£ 


pH H -x 
o 


+ 1) J dZ J dZ 


( 30 ) 


where 


£ - (Vl [« 3 - X) 2 - c 4 2 ] 


IV x I Sc 4 


f = 0 


X _ X > c 

3 I c 4 


(b) The case for rolling and sliding 

The rolling speed of the lower cylinder u ^ can be expressed 
in terms of u^ with an appropriate coefficient which depends 
upon a sliding speed desired between the two cylinders. 

Let = b^u^, and substituting b^u^ in Eq. (22) 

Bp ,;'48p)u N i . . P V i r° r 

3X 


H H“ 
o 


= cfs) t (0 - 5)(1 + b i> vv fi-irL L (0o5)(1 + v (f i> 


PH -X 
o 


a(PH o ) . - 
BT 


+ f + (1 - v< x 3- x) \-=r-) p j dx ; 


8P 


HZ 




(31) 


Q = \ 


/'48Ua 


) r U°- 5 > (^r 1 ) ( H r rv - tv) f i (°- 5 

H 


2 ‘ 

H -°° 
o 


pH H -X 
o 


( 1 + V (lf: + VoV^ 


+ f — + 
BT ST 


(i - b x )(x 3 - z) (-55-) pj dz} dx 


(32) 


where 
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f(x,x 3 ) = (^-) j_ (X3 — x) 2 - c 4 2 j 

r 


f(x,x 3 ) - 0 


The dimensionless load becomes 


1- r*° 


W = 4P 


HZ L 


P(X jT)dx + 


n] 


|x 3 - x | s c 4 

|x 3 - x| > c 4 


(33) 


2.6 Method of Solution 

2,6,1 Outline of Approach 

Since the pressure profile in the second half of the contact 
region is assumed to be the Hertzian pressure profile, the solu- 
tions for the first half pressure and film thickness profiles are 
obtained in the present study. The first half of the contact 
region is divided into two regions - the inlet and middle region. 

The increase in the inlet pressure is gradual and does not 
reach a high value comparing with the pressure in the middle 
region. Thus, the nonlinearity of the governing equation is not 
severe, and consequently the direct iteration method pan be used 
for the calculation of inlet pressure without introducing any 
convergence difficulty. Rather than calculating pressure directly 
from the governing equation, Q is obtained first by Eq. (30) for 
pure rolling case or Eq. (32) for rolling and sliding case. Then, 
from Q the inlet pressure is obtained. 

In the middle region, the system equations in the finite dif- 
ference form are solved by the Newton-Raphson method. The solution 
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of the system equations gives the pressure correction terms at 
every grid point in the middle region. Since finer grid spacings 
are required near the asperity to account accurately the variation 
of pressure and film thickness around the asperity, it is necessary 
to change the grid spacings for successive time intervals as the 
asperity moves toward the center. Thus, at each time interval, 
the previous values of pressure, film thickness and density are 
determined by linear interpolation to use in the calculation of 
time derivatives of these variables. Details of numerical treat- 
ment for the pressure and film thickness profiles and the center 
film thickness will be given in the next section. 


2.6.2 Integration of Pressure in the Inlet Region 

The equation used for inlet pressure calculation is written 

at Kth grid point and time T • 

m 


.... - w_ f{ ty (■, ■ h - ty r t- 


pc r S(pH o ) 


m 


r m 


pH H m -X 
o 


ST 


+ f (i + i) J dz ) - <V <V 

m 


(30) 


where 


U 1 ^ 


48Ua) 

V 

H ™ 
o 


In evaluating the integral in Eq. (30), the pressure, density and 
deformation are considered to be known and are taken as the pre- 
viously ihterated values. 
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Defining -x^ as th e dividing point between the inlet and 


middle region, Q . can be written as 

KAjin 


Q KA,m = (U 1 )(<Z W 


From the above equation 


T CAjm 

1 «kA 


(34) 


Substituting Eq. (34) for into Eq. (30), one can eliminate 
in Eq. (30) 

(35) 

t 

where Q is determined by the solution of the system equation 
KA,m 

in the middle region. Once Q is determined, through Eq. (15) 

ix y in 

P T . is obtained. 

K,Itl 

For the case of rolling and sliding, Eq. (32) is used to 
obtain the inlet pressure following the same method described 
above. 


2.6,3 Calculation of Deformation 

The numerical quadrature for the singular kernel in Eq. (30) 
is the same as the one detailed in Section 2.5.3 in Part I. How- 
ever, the upper limit of the deformation integral in Part II 
differs from that of Part I. The pressure distribution in the 
second half of the contact region is the Hertzian profile and the 
width of it is also the Hertzian half -width a. 
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1 Kf-2 

I V’ 4 "! z - x K l dz -I [ p j,„ K i(- x K,- z j ) 

'Hi j=i ,3,5 


+ P.., K 0 (- X - Z. ) + P K~ ( - X - Z . J 1 (36) 

j+I,m 2' K, j/ j+2,!n 3\ K, j/ J 


1 Kf-2 

I P ,n<« 4 "l Z l dZ ■ 1 [ P J ,„ K l (- X K0,- z j) + Vl.nH X K0,- V 

KI j=l,3,5 


+ P j+2,m K 3(' X K0,’ Z j) ] 


(37) 


where K^, ^ and were defined by Eq. (42), (43) and (44) of 

Ref. 6, and X is unity. When Z. is larger than Z - Z . is 
K.r j KO , j 

replaced by Z^ in Eq. (42), (43) and (44) of Ref. 6. 

In order to facilitate the differentiation P„ with respect 

K,m r 

to P . K, , K. and are rearranged in such a way that P. has 
j,m V 2 3 6 J j,m 


a single coefficient R(-X v - Z.): 

J 

R (- - S l( - Zj ) 


j - 1 


■ v- v Z j- 1 > 


even 2 ^ j 5 Kf-1 


= S 


3 ( ' V Z j-2 } + S l ( ' V' V 0dd3 ^ 3 S Kf-2 


S 3<- V Z i-2 } 


j = Kf (38) 


where 


S (- X Z.) = K (- X Z.) - K (X - Z.) 
n K, j n K, j n KO, j 


(39) 
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The final form of the deformation equation is 


K,m 


KO 

- C, y R(- X„ - Z.)P. 

5 ^ K> J J,m 

j=l,2. 


where 


c -ifle. 

5 H TT 


( 40 ) 


(41) 


2.6.4 Elastohydrodynamic Equation in the Middle Region 
(a) The case for pure rolling 

The governing equation written at Kth grid point and time 


T 

m 


is 





„ / *» L 9T + 1 V3x + dT' J X J 

H m 


( 

'C 

p H 

K,m o 


m 


(42) 


K,m om 


(43) 


where 


H. = 1 + 
K,m 


16PH N , _ 

ir^-J X K + 

om 


K,m 


(44) 


_ r I6p Hz \ r (. 


■ „ 2'; 

"K,m V - J L V X 3 " X K > ' C 4 J 

r m 


|X 3 - X K | ^ C 4 (45) 


m 


f K,m = ° 


x„ - x > c, 

'3 K. 1 4 

m 
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The integrand in Eq. (42) can be split into two terms - pres- 


sure dependent and independent terms. 

. « ftp . 


^ o ; - /do do'' r - - 1 

ST + f K,m W + ST ) = ' L^ PI V + f K,m P K,m-lJ /AT rr 

K,m K,m-1 ’ 


r _ _ _ ”i 

+ | (PH ) + f_ 

1 o „ K,m K,m 


P„ | /AT + f 

„ K,mJ m k,„. ^ ,, 

K,m * K,m 


\sx J 


(46) 


K,m \dX. 

The first term in the right hand side of Eq. (46) is independent 

of P. and is defined as Y (- x„) . The next two terms are de- 
j,m m TC' 

pendent of P m and defined as T^(- x^) . The absence of deforma- 
tion term in the integral makes it simpler in obtaining the pres- 
sure derivative of the integral and removes the convergence dif- 
ficulty between two successive iterations. 

Using the trapezoidal rule, the integral (42) can be written 
as : 


KO 


r° r 3(P H 0 ) _ 'i i “ r i 

J L— 55— + £ ^ * fii J d * * 5 L LV V + 1<- V j 

" i-K 


= I, 


K,m 


(47) 


where 


ffi i = X i+1 “ X i-1 


K+l^i^KO-1 


= X. - X. 
l+l i 


i = K, KO 


Eq. (30) is expresses in the discretized form at anc * 

T , where the pressure dependent variables are replaced by their 
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respective functions, 


T m ( P K+l/2/ ' PK,m ) l 1+<: 6 ) 4 +l /2- C 5 i R (‘ W‘ Z j) P j ,, 


2 

+ C 7 [ ( X 3 " X K+l/2 7 C 4 - / ( ex P^ a p ; 


K+l/2 ,m' 


J, 

[ ltc 6W-si R (- W,- z j) p j,, 


K+l/2,m 


I K+l/2,m 1 

/ B P„, . /<> \ / 

h u+ T , — ) 

O’”' 1 + A l P K+l/2,m^ 


16P HZ _ 8P HZ 48U 

H ’ C 7 - 1 and C 8 2 

om r H 

om 


Eq. (48) is one of the system equations written at an< ^ 

Similarily, we can obtain the system equations for the case 


of rolling and sliding. 


\ ( P K+l/2 ) ■ fc‘fA * ! l 1+C 64l/2-S I R V- Vi/2,- z j) P j,, 

* j-1 


“f 

V f 


2 .3 

+ C 7 - V X 3 " X K+l/2' ' C 4 J I " + b P ( ex P( a P k+ 1/2 ,m^ J 


,N T 

{ L 1 + °6 X K+l/2 ‘ C 5 2. R V" X ] 


: . z ^ p 

K+l/2 , j' j,m 


1 + A l P K+l/2,m 
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K+l/2,m 


( B P 

H o ,m \ + 1 + A 


K+l/2 ,m 
A l P K+l/2 , 


\ 

•J 


; 


(49) 

cont. 


Applying the Newton-Raphson technique to the system equations, 
we obtain 

{U p >} (n) + {ipj (n+1) <n) - 0 < 5 °> 

where j - and J represent a column matrix and a N x N matrix 

respectively, and A* indicates partial derivative is to be taken 

with respect to P . n is the level of iteration, 
m 

From Eq. (50) one obtains 

rl. 


t AP m} <,,+1) ' - 0'V P >] <n> - {U P >} <n> 


(51) 


The right hand side of Eq. (51) is assumed to be known from the 


lower level iteration, and 


~Fi 


(n+l) , 


is defined as 


W <n+1) - U <n+I) - k (n) 


(52) 


The elements of matrices in Eq. (51) are detailed in Appendix 


2.6.5 Calculation of Center Film Thickness 


From Eq. (30), the integrated Q T „ is obtained: 

Ko , m 
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48Ua 


xo ,m 


. {ftf (v^)-^ 

H m - “ H pm pH 

m o 


-) ■ 

3 J j 

pH H m -x 
o 


-d(pH o ) 

. 5t 



(53) 


Let 


r° i / p 

J — W - zr) dx = Q] 


-CO H 


m 


p m 


and 


rO /I \ 

J tsj * - q 2 

pH m 


Then Eq. (53) can be written as 


Q v H +(48Ua) Q.H +(48Ua) Q. = 0 
Xo,m om X 1 om ' x 2 


(54) 


Eq. (54) is the cubic equation of H . Since the coefficient in 

om 

Eq. (54) implicitily depends upon H Qm and P^, the analytical solu- 
tion for H Qm is impossible to obtain. Rather, using the so-called 

secant method, H is calculated numerically. 
o,m 


2.6.6 Outline of Numerical Procedure 

It is assumed that the center pressure is constant while the 
load and center film thickness change as the asperity moves toward 
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the center of the contact region. Except for the pressure distri- 
bution around the asperity, the pressure profile can be approxi- 
mated to the Hertzian profile of the same center pressure, which 
is used as a initial estimated value for pressure at the first 
time step. The calculation starts with the asperity located far 
away from the inlet of the contact region. 

Written below are the procedures of numerical calculation at 
each time step: 

1) At the first time step, the solution for the steady state 
EHD problem is used. This solution can be taken as a true 
solution for the transient EHD problem because the asperity 
located far outside of the inlet of the contact region can 
not influence the pressure distribution inside the contact 

region. From the second time st e p on, the initial es timated 

pressure is obtained by interpolating the previous pressure 
distribution for new grid spacings. 

2) Using the initially estimated pressure distribution, the 
film thickness, density and viscosity of lubricants are cal- 
culated. Also the center film thickness is determined using 
the pressure and film thickness profiles at the previous times 
but the new asperity location in the calculation is incorpo- 
rated. Then the system equations for the calculation of pres- 
sure correction in the middle region are solved by the Newton- 
Raphson technique. The inlet pressure is obtained by the 

linear-interpolation with the factor P,, A ^ where 

r KA,m KA,m 



P„. is obtained from the system equations, and then the 
KA,m 

film thickness is calculated using the newly obtained pressure. 
3) If the converged solution for the pressure in the middle 

region is obtained, the inlet pressure is recalculated by 
Eq. (35) and at the same time the center film thickness is 
also determined. At this time the overall pressure distribu- 
tion is checked for convergence. If it is converged, the 
load is determined by Eq. (33). Otherwise, the proce- 
dures (2) and (3) are repeated until the converged solution 
is obtained. 

The computer flow diagram and listing are shown in Appendix C. 
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CHAPTER 3 


DISCUSSION OF RESULTS 


3.1 Introduction 

Since the present study is mainly concerned with the variation of 
level of the separation between two cylinders, the results are presented 
as series of film thickness profiles with corresponding pressure dis- 
tributions as the asperity moves toward the center of the con- 
tact region. The asperity height is varied from 1/2 H g to 2 h^ and the 
asperity width is varied from a/4 to a/2, where H g is the steady 

state, center film thickness and a is the Hertzian half-width. The 

-12 5 

conditions used in the present study are: U = 5.3 x 10 , p q = 10 psi 

and a = 0.95 . 

f 



3.2 Pressure Profile 


The steady state pressure profile and the corresponding film thick- 
ness profile in the absence of the asperity is shown in Fig. 2. These 
profiles are used as a reference to compare with the transient solutions 
obtained with the asperity under the same conditions. The pressure 
profiles in Fig s , 3 to 11 show the change in pressure caused by the 
asperity as it moves from far outside of the contact region to the center. 
In each series, only three pressure profiles are presented, other inter- 
mediate profiles calculated at time intervals between the three positions 
have been omitted. 

In general, the film pressure around the asperity tends to increase, 
and the increase becomes the largest when the asperity enters the con- 
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tact region. The increase in pressure is very closely related to the 
squeezing action of the asperity. The magnitude of the squeezing 
action increases as it approaches the contact region, and reaching a 
maximum at the inlet of the contact region. It then starts to 
diminish as the asperity moves further toward the center. It is evident 
that the greater the squeezing action, the greater the pressure increase. 
The pressure fluctuation seems to be insufficient to cause any adverse 
effect on the contact surface. The elastic depression on the base sur- 
face of the asperity did not occur, which is in the direct contrast of 
the results [ 7 ] in which the pocket is formed elastically when two one- 
dimensional asperities approach each other. 

When the asperity approaches the inlet of the contact region, the 
inlet pressure gradient becomes very steep. The inlet pressure increases 
from the ambient value to a very high pressure in a short distance. The 
lubricant behind the asperity is less pressurized while the lubricant 
in the front of the asperity is severely pressurized by the squeezing 
action of the asperity. 

Shown in Figs. 10 and 11 are the results of the rolling and 
sliding of the two cylinders. Fig. 10 is the result for = 0.9U 
and = 1.1U, and Fig. 11 is the result for = 1.1U and = 0.9U. 
Both pressure profiles have the small local pressure peak when the 
asperity is in the contact region. As mentioned previously, 
the rolling term and the squeezing term of the asperity do not cancel 
out each other. Thus, even when the asperity moves in parallel with 
the flat surface, the squeezing action is still possible allowing 
positive and negative pressure on the lubricant in either side 
of the asperity. However, when the asperity is near the inlet of 
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the contact region, the squeezing action overcomes the sliding term 
and consequently the pressure bump is not generated in the fluid film. 
When the asperity moves faster than the lower cylinder, the pressure 
in the left side of the asperity changes more rapidly, and the reduc- 
tion in pressure is quite substantial compared with the pressure pro- 
file with the asperity entering the inlet of the contact region. 

3.3 Film Thickness 

The steady state film thickness profile in the contact region is 
approximately parallel with the flat surface, and the thickness is de- 
pendent upon a rolling speed if other conditions are the same. Since 
the conditions in both steady and transient problems are the same, the 
quantitative comparison between them can be made with regard to the 
effect of an asperity. 

— — — The- shape-of-the--transient _ film _ thickness prof Tie - is notably HTf- _ 
ferent from the steady state profile. The steady state case profile is 
approximately constant across the contact width, whereas the transient 
film profile is sloped considerably toward the contact center. The 
slope of the film thickness profile increases with increasing asperity 
hieght. This phenomena helps in preventing the direct contact of the 
asperity tip and the flat surface. As far as the shape of the film 
thickness profile is concerned, the influence of width of the asperity 
is not significant even though the wider asperity tends to increase the 
level of the contact separation significantly. 

In all cases studied, the center film thickness is always larger 
than that for the steady state case. The increase in center film thick- 
ness is dependent upon the asperity height, and the level of separation 
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increases with increasing height of the asperity. Even for the asperity 
with the height equal to twice that of the steady state film thickness 
H g , the tip of the asperity does not touch the flat surface when it passes 
through the first half of the contact region. This means that the level 
of separation increases more than 100% of H . This result may be one 
of the most significant findings in the present study and appears to 
lend credence to previous speculations regarding the beneficial effect 
of the asperity on the film thickness. 

Displayed in Fig. 12 is the center film thickness vs. the asperity 
location curves, which shows that the center film thickness increases 
continuously as the asperity moves toward the contact center. The rate 
of increase in the center film thickness is dependent upon the width 
and height of the asperity. For the same height of the asperity, the 
center film thickness for a wider asperity increases faster and larger 
than that for a narrower asperity. When there is sliding between the 
two cylinders, the change in center film thickness depends upon a speed 
of the cylinder to which the asperity is attached. In the present study 
the asperity is attached to the upper cylinder. The results show that 
if the upper cylinder moves faster than the lower cylinder with the same 
rolling speed of pure rolling, the center film thickness increases more 
than the one for the case of pure rolling having the same asperity geo- 
metry. In the opposite case, that is, when the lower cylinder moves 
faster than the upper cylinder, the increase in center film thickness 
is less than that for the case of pure rolling. The faster the upper 
cylinder can pressurize the lubricant more effectively in the contact 
region. Consequently, the level of separation tends to increase to ac- 
commodate the lubricant swept in by the asperity. 
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Fig. 13 shows the relationship between the inlet minimum film thick 
ness and the asperity height, both of which are normalized by the steady 
state film thickness to show the extent of the increase of the center 
film' thickness by the asperity. The dash line in Fig. 13 indicates the 
relation between these two values if the center film thickness is not 
raised by the hydrodynamic effect of the asperity. For small asperity 
heights (f/h g <<: 1 ) , the solid lines coincide with the dash line indicat 
ing an absence of hydrodynamic effect due to the asperity. For large 
asperity heights (f/h s > 1) the dash line indicates a negative minimum 
film or an interference, whereas the solid lines still show a clearance 
between the tip of asperity and the opposing surface. Even for the 
very severe case of f/h g = 2, a clearance of approximately 10-20% of 
the steady state film thickness was found to exist underneath the tip 
of the asperity. 


3.4 Load 

It was found that the load carrying capacity of two cylinders with 
an asperity is larger than that for two cylinders without an asperity. 
However, the increase is not more than 15% of the steady state load with 
the same peak pressure. Furthermore, in the load calculation the pres- 
sure distribution in the second half of the contact region is assumed 
to be a Hertzian pressure profile. Thus, it may be possible that the 
load may change unfavorably when the asperity enters the second half 
of the contact region, which remains to be investigated. 
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CHAPTER 4 


SUMMARY OF RESULTS 


The important aspects of elastohydrodynamic lubrication, with a 
single, one-dimensional asperity , have been found by solving numerically 
the coupled transient Reynolds equation and the elasticity equation. 

Even though the assumption of a single asperity is highly ideal, but 
this study sheds some light on the effect of surface roughness on 
elastohydrodynamic lubrication. 

The results show that: 

1) The film pressure tends to increase more than the 
steady state pressure, and in particular, the increase in 
pressure reaches a maximum as the asperity approaches the 
inlet of the contact region. The asperity height and the 
pressure increase above the steady state pressure are closely 
related to each other; the higher the asperity height, the 
larger the pressure increase. In the pure rolling case, it has 
been found that a local pressure peak is not developed. How- 
ever, in the cases of sliding and rolling, a small, local 
pressure peak is developed on the pressure profile when the 
asperity moves into the contact region. 

2) In general, the overall film thickness profile increases 
with increasing asperity height, but is not significantly af- 
fected by the asperity width. Moreover, the slope of the over- 
all film thickness profile for the transient cases is much 
greater than the steady state profile which is approximately 
constant across the contact width.. The increase in the center 
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film thickness also depends upon the width and height of the 
asperity. Even for the case of an asperity height of 2H g , 
the center film thickness increases more than 100% compared 
to the steady state center film thickness. 

As mentioned before, the surface condition employed in the present 
study is highly ideal. Thus, the present results may not be applicable 
to a more realistic surface condition of randomly distributed asperities. 
However, the results of the present study suggest that the rough contact 
surface is beneficial in generating continuous fluid film between the 
heavily loaded, two contact surfaces. 
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Straight Exponential Lubricant Model 



Fig. 2 Pressure and film thickness profiles without an asperity, 
pure rolling. 
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Fig. 3 Pressure and film thickness profiles with an asperity 
pure rolling, asperity width = a/4 and Height = H /2. 






Fig 


5 Pressure and film thickness profiles with an asperity 
pure rolling, asperity width = a/4 and height = 1.5H 



Fig. 6 Pressure and film thickness profiles with an asperity, 
pure rolling, asperity width = a/4 and height = 2H g . 
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Fig 


7 Pressure and film thickness profiles with an asperity 
pure rolling, asperity width = a/2 and height = H /2. 



Fig. 8 Pressure and film thickness profiles with an asperity, 
pure rolling, asperity width = a/2 and height = H g . 
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Fig. 10 Pressure and film thickness profiles with an asperity, 

rolling and sliding, asperity width = a/4 and height = H , 
U 1 = 0.9U, U 2 = 1.1U. 
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Fig. 11 Pressure and film thickness profiles with an asperity 
rolling and sliding asperity width = a/4 and height = 
IL = 1.1U, U 9 = 0.9U. 


Asperity Width = a/4 With No Sliding 

. — Asperity Width = a/2 With No Sliding 

■ Asperity Width = a/4 With Sliding 
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12 Center film thickness vs. asperity locations. 
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APPENDIX A 
SYMBOLS LIST 

Half of Hertzian width 

coefficient of density 

coefficient of density 

constant in deformation formula 

constant in deformation formula of cylinder 1 

constant in deformation formula of cylinder 2 

coefficient of deformation formula 

half width of asperity 


Deformation 



Equivalent Young' modulus 
Young's modulus of cylinder 1 

Young's modulus of cylinder 2 


Height of Asperity 


G = OE 


h 


h 


1 


Film thickness 
See Eq. (7) 
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Rigid center film thickness 



center film thickness 
geometrical film thickness 

Minimum film thickness 


A dummy index 
See Eq. (B.7) 

A dummy index 
A dummy index 
A index for time step 
Pressure 
Center pressure 

Hertzian pressure 
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r 


See Fig. 1 


— r 

r = i 


R = 


R 1 x R 2 
R x + R 2 


“k,J 


T =- t 


S . = — j— (An|u. 
J 2 j 


U j = ' Z j' X K 


- U l- 


U = 


Vl 

ER 


Radius of equivalent cylinder 
Radius of cylinder 1 

Radius of cylinder 2 

See Eq. (A. 10) 

time 

For Part II 



Speed of upper cylinder 
Speed of lower cylinder 
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X 


X 


3 



X 


3 



X 


KA 


*KA" 




coordinate along film 
Instantaneous location of asperity 


Coordinate separating the inlet and middle region 


Load per unit width of cylinder 


Dummy coordinate along film 


Poisson's ratio of cylinder 1 
Poisson’s ratio of cylinder 2 
Density 

Ambient density 


viscosity 
Ambient viscosity 
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a Pressure-viscosity coefficient 

a = o;P 

o 


B 

? = B 


Second pressure-viscosity coefficient 


V " 0 


A-Y (p) 
m 




System equation 

Derivative of Y (p) with respect to p 
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APPENDIX B 


CALCULATION OF MATRIX ELEMENTS IN EQ. (50) 


[&.Y (P)] in Eq. (50) has N x N elements, each one of which is the 

derivative of Y (P) with respect to P . For convenience. Y (p) is 
m * m m 


written below. 


T „<W> - (• 


P - P 

K4-1 ,m K,nT 


*=) {l + C 6 4 +1/2 - C 5 y R (- X k+1/2> - Z ) P jjt 

T J-l 

+ c 7 [ (x 3 - X R+1/2 ) - C 4 J } - (Cg) (exp (a P K+1 / 2jin )) 


{[ 


K f 

1 + C,x:., d R l- X, 


J 6 K+l/2 5 


K+l/2 


J-l 


- Z. ) P. 

» 3' I,™ 


K+l/2, m 


l + ^>±l/2, m . J H ( 1+ ,- B W» m ) 


(48) 


1 + A P 

A 1 K+l/2 ,m 


om y 1 + A,P . ,, 

1 K+l/2 ,m 


Eq. (48) is EHD equation in which viscosity, density and film 
thickness are expressed as a function of pressure. Before differentia- 
tion, the algebraic average of these variables is identified at - X^^^ > 
and Eq. (48) is expressed in the following form: 


= <I> (\+l,t 


vw> 





(h + 

K+l ,m 


K,m 



P K+l/2 




,m 


(A-l) 
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- «S> ) (i K+1 ,„ + } 


H om P K+l/2 ,m 


CA-1) 

cont. 


where 


KO 

v»' ( t> l [V- V + V- VJ “i 

i=K 


and the film thickness m an£ * C ^ e ■*- nte 8 ra ^ P k+ 1/2 m are ex “ 

pressed as the average of the two values at -X^ and as: 

I K+1 / 2 , m = 2 ( I K+l,m + ^.m) (A_2) 

^+1/2, m “ 2 (nR+l.m + ^.m) (A " 3 ^ 

(A-4) 


H 


= \ (h + H ) 

A K+l/ 27 m --K+lTtn — K T m— 


Hr+ 1/2 m t * le thickness including the height of the asperity and 

H 1 is the film thickness excluding the height of the asperity. 

K+l/2,m 

p__ and pi _ are assumed to be a function of the average 


pressure, j (P K+1 _ n + as: 


1 ! \ 

2 B V P K+l,m f P K,m^ 

J 


K+l/2 ’ m l + k (p^ + P K 

2 1 \ K+l,m K,i 


(A-5) 


^K+l/2 ,m eXp l_2 “ V P K+l,m + P K,m/ -i 


(A-6) 


Eqs. (D-2) to (D-6) are differentiated with respect to P m> where 

H_, and H, are the functions of P. regardless of indice j 

K+l/2,m 1 J,m 

K+l 3 tn 
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and the rest of Eqs. (D-2), (D-5) and (D-6) are the function of P. 

J 

only for j = K or K+l. The derivatives of these equations are written 
as : 


5H K+L/2,m _ 


BP. 
J 


^“k.j* V- +H k,J 


(A-7) 


where 


RRj, , = R(- \ +1) ~ z j) + R(- \ - Z ) for j = KA, KA+1 

*Yi ^ 


KA 

v { [ R (- W h ) + R< - V- V-i i 




for j = KA, KA+1 


To account for the effect of inlet 

pressure distribution on D„, + D.,,,. , the sum of the product of the 

v KA,m KA+1 ,m 

deformation kernel and inlet pressure ratio is considered as the deri- 
vative of the deformations at -X^ and with respect to P^ ^ 


and P 


KA+1 ,m" 

For j = K, K+l , 


B|l, 


K+l/2 ,m 1 

BP. ' 2 ^K+l/2 m 

J 


( A-8) 


Bp K+l/2,m _ 
BP . 


B 

—r 


j,m 2 + A. \P . + P ) 

J K+l ,m K,m/ 


3 L, -i /H + f f \ /Bp \ 

- /.1a ( om j j t ^ j /AX ) 

■ V \ AT AX j \3p. / '•“V 

m K j ,m 


BP. '2 

j,m 


(A-9) 


(A- 10) 


For j / K, K+l, Eqs. (D-8), (D-9) and (D-10) is zero. 
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Using Eqs. (D-7), (D-8) , (D-9) and (D-10), the derivative of 


\ (P K+l/2 ) is written as: 


/ \ 


3 Y m^K+l/ 2 J P K,m^ , , e „ „„ s h .77 ^ 

“ST v ^ J ( - 1 * 5C 5 RR K,3 ) V\+i jm + 

i ,m K 


■« s g ><i> - <C 8 > (““*<«*: 


K+l /2 


{ [i (i _+I, ) - ^ j - I (_i ) ( W+ h j } 


K+l,m R,m p 


K+l/2 ,m 


^om^K+ 1 / 2 ,m 


■<c 8 ) (“P<“ Vl/2,.>) {<-°- 5C 5 E1, K, J > + 


P K+l/2,m P + *1 ^ P K+l,m + P K,m) 


2H om ^+1/2, m 2 + \ V P K+l,m + P K,m7 


~ 0 ( I K+l,m + ’k.J 


( 1 T ( /H om + ^K,m , ^K,m\ ( B 

- (, — = ; L > at — + b?~J 7 

2H p__ . . n m 2 + A. U> + P„ r 

om K+l/2, m 1 ' K+l,m K,r 


,) <«*> 


/H + 


om K-H,m , K+ljin’N f 


. / om K+i 3 m , 
+ V AT + 


gr 1 - 2 ) ( B ,) )] } (D - IX) 


where 


6 = 1 


j = K+l 


6 = - 1 


j = K 


Eq. (D-ll) is one of N x N matrix elements, the expanded form of 
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F.q. (50) is 




s w 

3Y (P ) 
m'' Klv 

3Y (p ) 
KA' 

3P KA,m 

SP KA+l,m 

5P KO-l,m 

a VW 

*VW 

a vw 

^ P KA,m 

1 

1 

| 

3p ’ 

KA+1 ,m 

1 

i 

Sp 

KO-l.m 

1 

1 

1 

1 

It 

a V P KO-l> 

1 

1 

VkO-P 

1 

1 

BY (P ) 

9P KA,m 

?P KA+l,m 

aP KO-l,m 


S. ^ 


/ x 

Ap 

KA,m 


AP„. , 
KA-1-1 ,m 


AP 

\ 


KO-l,mJ 




Y (p ) 

V KA ; 


y (p ) 

m v KA+]/ 


y (p ) 
V KO-l ; 


/ 


The inversion of the above square matrix is obtained by the Gaussian 

elimination method. {AP t , } is the column matrix and added to {p„ } 

K,m h,m 

at each iteration. When [A } does not meet, the convergence criteria, 
the iteration is repeated with the adjust inlet pressure by linear in- 
terpolation and a constant center pressure until the converged solution 
is obtained. 
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APPENDIX C 


COMPUTER PROGRAM FLOW DIAGRAM AND FORTRAN LISTINGS 



Fig. B-l. Flow Diagram For Program 
NASA 2. 
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n n 


PROGRAM NASA2 ( INPUT, OUTPUT .PUNCH, TAPE3»TAPE5=INPUT,TAPE6=OUTPUT ) 
COMMON K.I ,KO,K.B»KBF,KF,Ml ,MF,K.H,KK.I »X2 »KXO * KK.B »K.02 * K.KF » K.F2 » 
1K3»MMI .KKH.KBI , KBFF » I S,N * NN , XH , U2 ,PO » C2 » RC » X3 » PC * ALPH , E » P I I ♦ 
2P1»H.HH,P,DEN»DEND»VIS,VISD,DE1,H1 ,PAI ,vd,psi »psd,ps»denx,c, 
3:<W»I »J,K»XC,KB2,KB1,EN,ED,PH,M,EE,SUBV»KIB,KI6F,C1 »DX»S» IIS* 
4XA»Sl»HO»DT »U1 *R»Q»D»DP, DENT #C9»BET » D1 »NCT *C3»C4 *U »DEO * C6 »SOE ♦ 
5PT »PP »S3 »S4 

DIMENSION XH<65) »Rl AO. 40) »EE ( AO) »H(57) »H1 ( 57) »HH( 57) ,P( 65) , 
1P1(57) ,DEN< 57) ,DEND( 57) ,SOE< 57 ) , PAI ( 57 ) ♦ V I S ( 57 ) , V I SD ( 57 ) . 

2D (57) »PS( 57) »DX( 57) ,DP< 57) »VD( 57 ) * A ( 57 ) *S( 57) .Q( 57.65 ) * 

3DEK57J »H0( 21) ,SUBV(65 ),F(65),C(2.3»65)»BET(40»40)»D1(57)» 

4PP ( 57 ) » PT ( 57 ) 

DIMENSION TEM( 40 ) « NI(40#4) 

REWIND 3 

READ (5*100) KA. K0» KIB, XI BF » KF. MI » MF ♦ KH 

READ (5*101) H3 »P0 *C5 

READ ( 5 * 103 ) PE, EU, EP, PC 

READ ( 5 » 108 ) E, EN, ED, ALPH 

KB=KIB $ KBF =K I BF 

t>H = PO/E 

S4=l • 1 $ S3=0.9 
XO = 57 $ KF = 6 5 
G= ALPH* 336 • 0$ KI=1 
PI 1=3.141593 

U= ( H 3*0 .75 / ( 1.26*G**0.6*PH*» 1-0.27) ) )**( 10./7. ) 

HO ( 1 ) =H 3 

P ( K I ) =0 .0 $ DEN ( K I ) = 1 • 0 $ VIS(KI ) = 1 . 0 

C9*l. 0/8.0 $ C3=16.0*PH**2 S C4=C3/PII 

U1=48.0*U/(HO( 1 )**2) $ C1=C3/H0( 1 ) $ C2=C1/PII 

C6=C5/C9**2 

XC=C6/H3 

WRITE(6,11) PH,E,EN,ED,ALPH 
WRITE(6,12) U»G,H3»U1 ,C1 ,C2 
WRITE<6,13) PO ,C5 
HTM=HO( 1) 

XXI=KI+1 $ K2=K I +2 $ KK0=K0+1 $ XXB=KB+1 
K02=X0-2 $ KKF=KF-1 $ KF2=KF~2 S K3=KI+3 
XA1=XA-1 $ XXA=KA+1 
MM I = M I + 1 S KKH = K.H+1 
K01=K0-1 

QX0=( 1.0-1. 0/ EXP (ALPH) ) /ALPH 
XH(KI)=-4.0 $ DT=1. 0/16.0 
C DETERMINE GRID SPACING AT FIRST TIME STEP. 

DO 7 K=2 , 7 

7 XH(K)=XH(K-1 )+l. 0/4.0 
DO 6 K=8 , 13 
6 XH(K)=XH(X-1)+1. 0/8.0 
XH ( KO ) =0 ,0 

DETERMINE PRESSURE DISTRIBUTION IN THE SECOND HALF OF CONTACT RE 
-G ION BASED ON THE HERTZIAN PROFILE. 

DO 202 X=KXO, KF 
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202 XH(K)=XH!K-1 ) + l. 0/8.0 

00 201 K = K0» KF 

201 P!K)=SQRT! 1.0-!XH!K)**2) > 

WR I TE ( 6 ,210 ) ( K , P(K), K = KO, KF ) 

1 IS=2 $ IS=1 $ P ( KO ) = 1 . 0 $ K = KO $ CALL DHD $ DEO=DEN!KO) 
NCS=1 

NCT=1 

IF(NCT.EQ.I) 175. 176 

C READ DEFORMATION. DENSITY. PRESSURE OF PREVIOUS TIME STEP. 

175 READ! 5,623 ) ( D1 ( K ) ,K = KI , KO ) 

READ ! 5 » 623 ) (P1!K)»K=KI»K0) 

READ (5,623) ( DEI ( K ) , K=K I »K0) 

READ (5,623) ( PT ( K ) , K = K I , ICO ) 

READ ( 5 , 624 ) H0( 15 ) .Cl ,U1 »XC 
READ! 5,624) H0( 16 ) ,C1 ,U1 , XC 
C2=C1/P 1 1 

WR I TE ( 6 * 1 10 ) U1,H0( 1) ,C1,C2,XC 
WR I TE ( 6 ,668 ) !K»D1 ( K ) ,K = KI ,K0) 

WRITE! 6,668 ) l K » DE 1! K ) , K. = K I • K.0 ) 

WRITE! 6,668 ) ( K » PI ( K ) » K = KI »K0 ) 

DO 1714 1=1,16 

1714 READ ( 3 ) M, ( ( Q ( K» J ) » J = KI »KF ) »K = KI »K0 ) 

176 DO 200 M= 17 » 20 
LL=0 $ 10 = 0 $ LC = 0 
IF(M.LE.2) 561, 562 

DETERMINATION OF DIVIDING POINT BETWEEN THE INLET AND 
MIDDLE REGION, THIS POINT DEPENDS UPON THE LOCATION OF 
ASPERITY AT EACH TTME STEP'i — ~ 

561 KA = 39 $ 60 TO 570 

562 I F ( M.LE .4 ) 563, 564 

563 KA = 39 $ 60 TO 570 

564 I F ! M.EQ .5 ) 565, 566 

565 KA=37 $ 60 TO 570 

566 I F ( M.EQ *6 ) 567, 568 

567 KA = 35 $ 60 TO 570 

568 I F ( M. EQ. 7 ) 569, 571 

569 (CA = 3 1 $ 60 TO 570 - 

571 IFIM.EQ.8) 572, 573 

572 KA=33 $ 60 TO 570 

573 IF (M.EQ.9 ) 574, 575 

574 KA=28 S GO TO 570 

575 I F ( M. LE • 12 ) 576,577 

576 KA=27 $ GO TO 570 

577 K A=25 

570 IFIM.6E.18) 140, 141 

140 KO = 57 $ KF = 65 % KKF=KF-1 $ KK0=K0+1 $ KF2=KF-2 
<02 = KO- 2 $ K01 = K0-1 
DO 142 K=KO, KF 

142 p(K)=SQRT< 1.0-<XH( K)**2 ) ) 

141 CALL TRANS 
KBI=KB+4 $ KBFF=KB+20 

READ THE KERNEL OF DEFORMATION FORMULA WHICH WERE 
CALCULATED AND STORED IN MAGNETIC TAPE. 


57 



n n 


READ<3) M » (<Q{K.J)» J=KI. KF ) ♦ K=KI» KO > 

WR I TE < 6 » 100 ) M 

WR I TE ( 6 » 15 ) KB I , KBF 

WRITE (6. 60 3) ( < K » J .Q ( K , J > ♦ J=45 ,65 > ,K=45»57> 

WR I TE ( 6 .668 ) (K, XH(K), K = KI, KO) 

WR I TE ( 6 .666 ) < K » DX ( K ) » K=KI » KO ) 

WR I TE ( 6 » 100 ) KBI. KBF 
IF(M.EQ.l) 65. 66 

65 P ( KA ) =0 • 35 
P ( 35 ) =0 .08 

DO 61 K=35 . 41 

61 P(K)=P(KA)-(P(KA)~P<35) ) * ( XH ( KA ) -XH { K ) ) / ( XH ( KA > -XH ( 35 ) ) 
SCC=XH( 35)-XH( 1) 

DO 67 K = 1 ♦ 35 

67 p(K)-P(35)-l XH ( 35 ) -XH l K > )*P(35)/SCC 
KA1=KA— 1 
KCA=KO— KA+ 1 
DO 586 K=K I . KO 
586 PA I ( K 1 =0 *0 

DO 583 K=KI.KCA 
DO 583 J = K I . KCA 

583 LET ( K ♦ J ) =0 .0 
PIP=1.0-P( KA 1 ) 

DO 584 K=«KA,KO 

584 P(K) =P(KA1 )+PIP*SQRT( 1 . 0- ( XH ( K ) /XH { KA ) ) **2 ) 

CALL INTEG1(XH(KI ) .XH(KF) »2,P,KF.VALUE»IER) 

I S=2S I IS=2 $ CALL DHD 

30 TO 37 

66 DO 174 K=KI.KO 
DEN ( K ) =DE1 ( K ) 

174 H1(K)=1.0+C1*0.5*XH(K)**2+D1 (K) 

DO 156 K=K I , KO 

IF(K.GE.KBI.AND.K.LE.KBFF) 157, 158 

157 H ( K ) =H 1 ( K ) $ HH(K )=H1 (K)+XC*1 ( ABS(X3-XH(K) ) )**2-C9**2) 
30 TO 159 

158 H(K)=H1(K) $ HH!K)=H1(K) 

159 P ( K) =P 1 ( K ) 

156 CONTINUE 

273 WR I TE ( 6 .668 ) ( K , DE 1 < K ) , K = KI , KO ) 

WR I TE ( 6 .668 ) (K, Hl(K), K = KI » KO) 

rFIM.L E.2) GO TO 522 

HO( M) =HO(M-l )*2.0-HO(M-2 ) 

WR ITE( 6 ,126 ) HO ( M ) $ GO TO 523 

522 HO ( M ) =H0 ( M— 1 ) 

523 PK A= P ( KA ) 

CALL INTEG1(XH(KI ) ,XH(KF) .2.P.KF.VALUE.IER) 

C CALCULATION OF CENTER FILM THICKNESS. 

I F ( M.GT .3 ) 186, 187 
186 CMX= 1 . 0E-20 $ HDEL=HO( M ) *0.00005 
HMAX = HO ( M ) * ( 1 .0 + 0.5 ) 

HMIN=HO(M)*( 1.0-0. 5) 

MT = 1 5 $ HO I =HO { M ) 

NRONE SUBROUTINE DETERMINES CENTER FILM THICKNESS 
AT FIRST TIME STEP. 

CALL NRONE (HO(M) »CMX .MT.HDEL.HOI. HMAX »HM I N ) 
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U1=48.0*U/HO(M)**2 $ WRITE(6»110) Ul» HO(M) 

1 1 S= 1 $ CALL OHD $ CALL INTEG 
HOO=HO(M) $HTM=HO(M) 

IF(M.GE.3I 60 TO 187 
S(KI )=0.0 i K.A1 = KA-1 
DO 183 IC=IC I »KA 

18 3 7D( K ) = ( HI IC )— DEO/DEN < IC ) — PAI ( K ) / ( OEN ( iC ) *HO ( M ) ) )/HH(K.)**3 
DO 184 K=K I ♦ KA 

CALL INTEG2 ( XH ( IC ) »XH(IC+1 ) » 1 .VD • iCO .VALUE » I ER ) 

184 SUC+l >=SUC)+VALUE 
DO 185 K=K I »KA 
QtC=Ul*S(K) 

185 P ( K ) = -ALOG ( 1 • 0-ALPH*QIC ) / ALPH 
WRITE(6»210) (IC,P(K> .K = KI *K.A) 

187 IF(M.GE.7) 341,112 

341 PUNCH 623, ( D1 ( K ) ♦ K=K I , <0 ) 

PUNCH 623,(P1(K) » K=K I » ICO ) 

PUNCH 623, ( DEI (K) • K = K I • ICO ) 

PUNCH 623, (PT(IC) ,K = <I,KO) 

PUNCH 624, H0(M-1)«C1,U1*XC 
GO TO 342 

112 I F ( M*GE • 3 ) GO TO 109 
I F ( NCT • EQ« 2 ) GO TO 109 
1 1 S= 2 $ I S = 2 $ CALL DHD 
NC=(C0-KA+1 
DO 278 K=K I , ilC 

DO„2-18_J=.IC.L» NC 

278 BET( IC, J ) =0, 0 
DO 279 K=KA,KO 

279 PA I ( K ) =0.0 
PKA*P ( K A ) 

GO TO 264 

342 f’KA=P ( KA ) 

WRITE(6»210) ( IC, P ( IC ) , K = K I * ICO ) 

1 1 S=2 $ I S = 2 $ CALL DHD $ I I S= 1 $ CALL INTEG 
WR I TE ( 6 »642 ) ( K » PAI ( K ) * K = KI,<0) 

WR ITE ( 6 » 2 13 ) ( K*HH ( K. ) , K=KI, ICO) 

WRITE(6»211) (IC»H(K)»IC = ICl» ICO ) 

3KA= ( 1,0-1.0/EXP(P(ICA)*ALPH) ) / ALPH 

S(KI )=0.0 $ KA1 = ICA-1 

WR I TE ( 6 *642 I ( IC »PAI ( IC ) , K = ICA • ICO ) 

WR I TE ( 6 » 1 10 ) Ul.HOIM) ,C1 »C2 
DO 645 K=KK I fKOl 
645 PS ( K ) =P ( K ) 

GO TO 109 
37 DO 35 K = K I » ICO 
35 PS < K ) =P ( K ) 

PKA=P(KA) 

89 I F < IO+l.EQ.l ) GO TO 109 
144 MT = 15 $ HO I =HO ( M ) 

LCC=2 

IF(M.EQ.I) 161, 162 

CALCULATION OF CENTER FILM IHICICNES8 BAbED ON l HE 
CORRECTED PRESSURE DISTRIBUTION BY NEWTON-RAPHSON METHOD* 
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161 CMX=1.0E-18 $ HDEL=H0<M)*0. 00005 
HMAX=HO<M)*1.2 
HMIN=HO(M)*0.8 

GO TO 163 

162 CMX= 1 • OE-2 1 $ HDEL=HO(M)*0. 00005 
HM I N=HO ( M ) *0 • 5 

HMAX = HO ( M ) * 1 « 5 

163 CALL NRONE ( HO ( M ) » CMX »'MT , HDEL »HQ I *HMAX »HM IN ) 

HCM=HO(M) 

HO ( M ) = ( HTM+HO ( M ) )*0.5 
HTM=HCM 

U1=48«Q*U/ HO ( M ) **2 

Cl=C3/HO(M) $ C2=C1/PI I $ XC=C6/H0(M) 

1 1 S= 1 $ CALL DHO $ CALL INTEG 
109 VR I TE ( 6 » 1 10 ) HO(M) »C1»C2»XC*U1 
173 WRITE(6*110) HO < M ) »HO < M- 1 > » DT » C3 
PKA=P ( KA ) 

4-6 1 I S = 2 $ 1 1 S = 2 $ CALL DHD 
302 IF(M.EQ.I) 264, 263 

263 1 1 5 = 2 s CALL INTEG 

264 I S=3$ CALL DHD 

C CALCULATION OF MATRIX ELEMENTS FO THE SYSTEM EQUATIONS. 

DO 90 K = KA, KOI 

n=k-ka+i 

SOE ( K ) = < H ( K ) +H ( K+l ) )*0.5 

A(K)=(HH(K+1 )+HH(K) ) *0. 5 

C8=( P(K+1 )-P(K) ) * ( A ( K ) **2 ) / ( V I S ( K ) *DX ( K ) ) 
IF(K.GE.KBl.AND.K.LE.KBFF) 411,412 

411 ST0=(S4-S3)*XC*( X3*XH(KBFF)-0.5*XH(KBFF)**2-X3*XH( K) 
X+0.5*XHJK)**2) 

GO TO 413 

412 STO=0.0 

413 EE(N)=C8*A(K)-U1*(0.5*(S4+S3)*(SOE(K)-DEO/DEN(K) )-STO 
X-0 . 5 /HO ( M ) * ( PA I (K+D+PAI (K> ) / DEN ( K ) ) 

224 DO 90 J=KA , KOI 
I I = J-KA+1 
QQ=0 .0 

521 IF(J.EQ.KA) 91, 92 

91 DO 93 I =K I * KA 

93 QQ=QQ+(Q(K, I )+Q( K+l , I ) } *P < I ) /P ( KA ) 

GO TO 99 

92 QQ=Q(K,J)+Q(K+1, J> 

99 R(N, I I)=-QQ*C2*(C8*1.5-u.25*li4+i3)*ul )+u.5»ul/Hui m)*iBE I iN+1 , 
XII ) +BET ( N , I I ) ) /DEN ( K ) 

508 IF(J.EQ.K) GO TO 94 

I F ( J.EQ.K+1 ) GO TO 95 
GO TO 90 

94 5 1 GN=-1 . 0 $ GO TO 98 

95 S I GN= 1 • 0 

98 RIN* II )=R(N» II )+C8*A(K)*(SIGN/(P(K+l)-P(K) ) -V I SD ( K ) / V I S ( K ) ) - 
X ( U1*DEND( K ) /DEN ( K ) **2 ) * ( DEO*0 . 5* ( S4+S3 ) +0 . 5 /HO ( M ) * ( PA I ( K+l ) 
X+PAI < K ) ) > 

90 CONTINUE 
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280 IG=K0~KA 

I E=K01-KA+ 1 

WR ITE { 6 » 1 10 ) C2, DEO. U1 

C SUBROUTINE INIS P IS THE OPERATION OF MATRIX INVERSION. 

CALL IN1SP(R.IG,1.E-7.IEER. 40. TEM.NI) 

IF(IEER) 153.153,154 
154 wR ITE ( 6 » 100 ) IEER 
GO TO 1000 
153 DO 105 N= K. I ♦ IG 
AS— 0 • 0 

DO 106 I I*K I » IE 
106 AS=AS+R(N,1 1)*EEU I) 

K=N+KA- 1 

DP(K)=-AS 

P(K)=P(K)+0P(K) 

PS ( K. ) =P ( K ) 

IF(ABS(DPUU 1-0.6) 105, 105. 503 
105 CONTINUE 

WR I TE ( 6 ,62 1 ) ( K. , DP ( K. > ,K = KA,K01 ) 

GO TO 504 

503 WR I TE ( 6 .62 1 ) (J , DP l J ) , J = KA »K) $ GO TO 1000 

504 WR I TE ( 6 ♦ 100 ) LL, 10 
PKA=P(KA) 

PW=0.0 $ PQ=0.0 
DO 171 K=KA » KOI 
PW = PW+DP ( K. ) 

PQ=PQ+P ( K. | 

171 CONTINUE 

PQ R sAB'STP W7P Q ) — 

KA1=KA— 1 
LL=LL+1 

128 IF(PQR.LE. 0.0005) 166, 48 
48 IF(LL.LE.IO) 301, 322 
301 DO 402 K=K I » K.A 1 

P(K)=P<K)*P(KA)/PKA 
402 PS(K)=P(K) 

547 PKA=P(KA) 

GO TO 461 

C INLET PRESSURE CALCULATION BY INTEGRATION . 

166 DO 241 K.=K I ,KO 

IF(K.GE.KBI.AND.K.LE.KBFF) 414, 415 

414 ST0=(S4-S3)*XC*(X3*XH(K.BFF)-Q.5*XH(KBFF )**2-X3*XH( K ) +0 . 5*XH ( K ) * 
X*2) 

GO TO 241 

415 ST0=0.0 

241 VD(IC) = ( (S3+S4)*0.5*(H< K ) -DEO/DEN ( K ) ) -STO-PA I ( K) / < HO ( M ) *DEN ( K. ) ) ) 
X/HH< K)**3 

QKA=(1.0-1.0/EXP(ALPH*P(KA) ) )/ALPH 
S ( K I )=0.0 $ KA1=KA-1 
DO 73 K=K I »KA 

73 S(K+1)=S(K)+0.5*<VD(K)+VD(K+1) )*DX(K) 

00 75 K=K I , (CA1 
P1(K)=QKA*S<K)/S<KA) 

QCC=ALPH*P1 (K) 

IF(QCC) 532, 533, 533 
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532 

533 
545 

75 

537 

534 
536 

322 

113 

115 

49 

116 

114 

652 

651 


WR I TE ( 6 *668 ) ( K , S ( K ) * K=K I » KA ) 

WR ITE ( 6 *668 ) ( K ♦ VD ( K ) » K = KI » KA ) 

LCC= 1 $ GO TO 534 

IF(QCC.GE.I.O) 545* 75 

WR I TE < 6 *66 8 ) ( K » VD < K ) . K = K I ,KO) 

KA=K-1 

WR ITE ( 6 » 100 ) ICA 
LCC= 1 $ GO TO 1000 
CONTINUE 
DO 537 K=K I ♦ KA 1 

P ( K ) =~ALOG ( 1.0-ALPH*P1(K) >/ALPH 
GO TO 322 

THE CALCULATION OF INLET PRESSURE BY THE LINEAR INTERPOLATION 
WHEN THE CONVERGED SOLUTION IS NOT OBTAINED FOR THE 
PRESSURE IN THE MIDDLE REGION. 

DO 536 K=K I » KA 1 
P(K)=P(K)*P(KA)/PKA 
PKA=P( KA ) 

PW=0.0 $ PQ = 0 • 0 

THE OVERALL CONVERGENCE TEST. 

DO 113 K=KK I ,K01 
PQ=PQ+P ( K ) -PS ( K ) 

PW=PW+P ( K ) 

LL=0 $ PQft=A3S ( PQ/PW ) 

IF(PQR.i-E. 0.0005) 651, 115 
10 = 10+1 

IF ( 10. LE. 10 ) 49. 651 
DO 116 K=KK I * KOI 
PS(K)*P(K) 

I I S=2 $ I S = 2 $ CALL DHD $ IIS=1 S CALL INTEG 
LC=0 

GO TO 144 

I F ( IO.LE.l ) GO TO 109 
DHC= (HOO-HO(M) )/HQ<M) 

IF(ABS(DHC) .LE. 0.002 ) 651* 652 
HOO=HO ( M ) $ GO TO 109 

CALL INTEG2 ( XH( Kl ) *XH ( KF ) » 2 * P , KF » VALUE ♦ I ER ) 


W= VALUE 

C1=C3/H0(M) $ C2=C1/PI I S IIS=2 $ IS=2 $ CALL DHD 
HTM=HO(M) 

WS=W*4.0*< PH**2 ) 

CALCULATION OF THE CENTER FILM THICKNESS BY THE 


NEW PRESSURE DISTRIBUTION. 

CALL NRONE ( HO ( M ) * CMX ,MT , HDEL * HO I *HMAX »HM IN ) 


U1=48.0XU/HO(M)**2 


WR I TE ( 6 * 220 ) 
WR ITE < 6 *210 ) 
WRITE16.213) 
WRITE(6»211 ) 
WRITE(6*215 ) 
WR I TE ( 6 * 126 ) 
WR I TE ( 6 *642 ) 
WR I TE ( 6 *643 ) 
DO 580 K=K I , 


M, 

( K » 

(K* 

(K, 

( K » 

WS 

(K.PAIIK), K=K I » KO ) 
C 1 » C2 ,H0( M ) 

KO 


W, HO ( M ) 

P ( K ) * K=KI, KO) 
HH ( K ) . K=K I » KO) 
H { K ) , K=K I » KO) 

D ( K ) » K = tCI » KO) 
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580 HH(K)=HH(K)#HO(M) 

WR ITE < 6 *21 3 ) ( X » HH ( X ) » X = XI » XO ) 

IF(M.EQ.I) 452. 453 

452 DO 454 X=XI ,X0 
454 PP t X ) =P ( X ) 

GO TO 456 

453 DO 451 X=X I »K0 
PS ! K ) =P ( X ) — P T ( X ) 

451 PP ( K ) =PT ( X ) 

WR I TE ( 6 .626 ) ( X»PS( X ) »X = XI * XO ) 

456 DO 182 X=XI » KOI 
SM=P(X)-1.05 
I F ( SM ) 182 » 182, 1000 
182 CONTINUE 

I F ( M« LE • 8 ) 200, 218 
218 WRITE(6,601) M,WS,HO(M) 

200 CONTINUE 
1000 STOP 

11 FORMAT ( 5H PH=,E14.6,3H E=,E14.6,4H EN=,E14.6,4H U1=»E14.6. 

X6H ALPH=,E14.6) 

12 FORMAT ( 5H U=,E14.6,3H G=,E14.6»4H H3=»E14.6»4H U1=»E14.6» 

X4H C1=,E14.6,4H C2=,El4.6> 

13 FORMAT ( 6H P0=,E14.6,4H C5=,E14.6) 

15 FORMAT ( 6H XBI=,I3,6H k8F=, 13) 

100 FORMAT (815) 

101 FORMAT(ElO.l.FlO.l.ElO.l) 

— 10 3-FORM AT <4E~10.-H ___ 

104 FORMAT ( F10, 1 , E10.1, ElO.l, E10.2, F5.2) 

108 FORMAT ( 3E10«2, F5.1) 

110 FORMAT ( 6E 1 5 « 8 ) 

126 FORMAT ( 6E15.6, 215) 

213 FORMAT! 3X//50X, *HH < X , M ) =*// 2X , *X* , 20X , *X* , 20X ,*X* , 
X2QX,*X*»20X.*X*,20X,*X*//<6( I3,E15.7,3X))) 

210 FORMAT ( 3X//50X, *P( K,M) =*//2X , *X* , 20X, *X*, 20X, *X*, 

X20X, *X*, 2 OX , *X*,20X, *X*//(6!I3, E15.7, 3X ) ) ) 

211 FORMAT (3X//50X, *H ! X«M ) =*/ /2 X » *X*, 20X, *X*, 20X, *X*, 

X20X, *X*» 20X, *X*,20X, *K*//(6(I3, E15.7, 3X ) ) ) 

214 FORMAT ( 10H SOE ( K ) = / / { 6 < I 4 , 2X , E 1 5 . 7 ) ) ) 

212 FORMAT ( 5H C8=,E12.5»7H EE ( X ) = » E 1 2 . 5 » 8H VI S ! X ) = . E 1 2. 5 
X , 8H DEN ( X ) = » E12» 5 ) 

226 FORMAT ( 5H CC=,E12.5.6H ROLL= , El 2 . 5 , 4H SQ=»E12.5.5H SUM=,E12.5) 

227 ( ORMAT ( 6H A ( X ) = , E 12 . 5 , 8H P A I ( X ) = . E 12 , 5 , 8H SOE ( X ) = , E 1 2 . 5 ) 

215 FORMAT! 3X//50X, *D < X , M ) =*/ t 2X , *X*, 20X, *X*» 2QX, *X* , 

X20X, *X*, 20X , *X*,20X, *X*//!6!I3» E15.7, 3X ) ) ) 

668 FORMAT (2X//(6! 14, IX, E 1 6 • 8 ) ) ) 

601 FORMAT ( 4H M=,I3,6H WS=.E15.7,7H HO l M ) = , E 1 5 • 7 ) 

603 FORMAT ( 3X//5QX, *MATR I X R ( X , J ) */ / ( 5 ( 2 1 4 , E14.5))) 

220 FORMAT ( 1H1 , 2X, *M=*, 13, 4X, *W(M)=*, E15.8, 4X, *HO(M)=*, E15*8> 
621 FORMAT! //3X»*DP< K ) =*// (6! I4,1X,E16,8>)) 

623 FORMAT ( 6E 1 2 • 5 ) 

624 FORMAT ( 4E20 • 10 ) 

626 FORMAT! 3X//40X,*THE PRESSURE DI FFERENCE*//2X,*X*,20X,*X* , 
X20X»*X*»20X»*X*»20X,*X*,20X,*X*//!6( I3»E15.7»3X>)) 

632 FORMAT !//3X,*EE(K )=*//! 6< I4,2X,E15.7))) 
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637 FORMAT (//3X**DEN(K)=*//(6( I4,2X.E15.7))) 

642 FORMAT(//3X**PAI (K)»*//<6< I4*1X*E15«7>>> 

643 FORMAT < 5H C1=,E15.6*5H C2=,E15.6»8H HO ( M ) = ,E 15.6 ) 

end 

SUBROUTINE TRANS 

THE CALCULATION OF NEW GRID SPACINGS AND DETERMINATION 
OF PRESSURE* DENSITY* DEFORMATION AT NEW GRID POINTS. 

COMMON K I ,KO*KB*KBF»KF,MI *MF*KH*KKI » K2 » KKO *KK.B * K02 » KKF » KF2 . 
1K3*MMI .KKH.KDI « KBFF * IS.N »NN ,XH *U2 ,PO » C2 . RC »X3 » PC ♦ ALPH » E , P I I * 

2P 1 *H »HH *P.DEN.DEND,VIS,VISD, DE 1.H1.PAI *VD»PSI . PSD , PS * DENX »C , 
3KW* I »J*K,XC*KB2.KB1.EN,ED*PH,M.EE,SUBV.KIB*KIBF»C1*DX,S» I IS* 
4K.A »S1 tHO.DT *U1*R»Q»D»DP. DENT .C9.BET *D1*NCT *C3»C4,U» DEO , C6 »SOE . 
5PT,PP.S3,S4 

DIMENSION XH(65)»R(40»40)*EE(40)*H(57)*H1(57) »HH ( 57 ) . P ( 65 ) * 

1P1 ( 57) , DEN I 57) ,DEND( 57 ) , SOE ( 57 ) *P A1 t 57 ) » V 1 S ( 57 ) . V I SD ( 57 ) ♦ 

2l>< 57 ) »PS(57 ) *DX( 57) »DP { 57) »VD(57 ) *A( 57) »S( 57) »Q( 57*65) » 

3DE1( 57 ) *HO( 21 > *SU6V(65 ),F(65)*C(2,3*65)»BET(40,40>»D1(57). 

4pP ( 5 7) » PT ( 57 ) 

MM = M 

KKK=0 $ K01=K0-1 
IF(MM.EQ.I) 55, 58 
58 L=MM/2 $ I F ( MM . EQ . L*2 ) 81* 82 

81 <B=KlB+M $ KBF=KIBF+M~1 $ KKK=1 
XBC=KB-2 

IF (MM.EQ.2 ) 71, 72 

71 XH< 14)=XH( 13 ) + l. 0/32.0 
XH ( KB )=XH( 14) +1.0/32.0 
GO TO 55 

72 DO 73 K=14, KBC 

73 XH(K)=XH(K-1 )+l. 0/16.0 
XH(KBC+1 )=XH( KBC ) +1.0/32.0 
XH(KB)=XH(KB-1 )+ 1.0/ 32.0 
GO TO 55 

82 KB=KIB+M-1 $ KBF=KIBF+M-1 * KKK=Q 
DO 74 K= 14 , KB 

74 XH(K)=XH(K-1)+1. 0/16.0 
GO TO 55 

55 KKB=KB+1 S KB20=KB+20 $ KKBF=KBF+ 1 
DO 12 K=KKB , KB20 

I F ( K • LE .KB+4 ) 16* 18 
16 XH(K)=XH(K-1 ) + l. 0/32.0 
GO TO 12 

18 XH(K)=XH(K-1 >+l. 0/64.0 
12 CONTINUE 

if(KKK.EQ.l) 24* 25 

24 XH ( KB+2 1 ) =XH ( KB + 20 > + 1.0/32.0 
XH ( KB+2 2 )=XH(KB+21)+1.0/32.0 
XH( KBF) =XH( KO+22 )+ 1.0/ 16.0 
GO TO 56 

25 KB21=KB+21 

DO 26 K=KB2 1 » KBF 

26 XH(K)=XH(K-1 )+l. 0/32.0 

56 CONTINUE 
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94 DO 96 K=KKBF » KO 

96 XH(K)=XH(K-1 ) + l. 0/16.0 
DO 97 K=KKO. KF 

97 XH(K)=XH(K.-1 ) + l. 0/8.0 
78 X3=XH( KB+12 ) 

DO 234 K=K I » KOI 
234 DX(K)=XH(K+1)-XH(K) 

113 kB1=KB-1 $ KB2=KB+2 S KBF1=KBF-1 

KB4=KB+4 $ KB16=KB+16 $KB17=KB+17 $ KB20=KB+20 
KB21=KB+2 1 

IF(NCT.EQ.1.AND.M.LE.17) 63.64 
64 IF(M.EQ.I) GO TO 63 
225 DO 67 K = KI. KO 

lF(K.LE.KB-2. OR. K. GE. KBF ) 68* 67 
68 DE1(K)=DEN(K)$P1<K)=P(K)$D1< K)=D(K)SPT(K)=PP(K) 
67 CONTINUE 

230 L=M/2 $ I F ( M.EQ. 2*L ) 201. 202 
201 DO 203 K=KB1 » KB2 

DE1(K)=DEN(K) $ P1(K)=P(K) 

PT ( K ) =PP ( K ) 

203 Dl ( K ) =D ( K ) 

DEl(KB+3)*DE,UKB+4) $ Dl (KB+3) =D< KB+4) 

PT ( KB+3 ) =PP ( KB+3 ) 

Pl(KB+3)=P<KB+3) 

DO 204 K=KB4 * KB16 
DE 1 ( K ) =DEN ( K+2 ) $ Dl t K ) =D ( K+2 ) 

PT ( K ) =PP ( K+2 ) 

204 Pi ( K ) =P ( K+2 ) 

DO 205 K=KB17 * KB20 

IF ( K.EQ.KB20 ) 219. 220" ~ ~ — 

220 L=K/2 $ I F ( K.EQ. 2*L ) 217. 218 

217 IFIK.EQ.KB17) 302. 303 

302 DE1(K)=DEN(K+1 > +0.5* ( DEN < K+2 ) -DEN ( K+l ) ) 

Dl (K)=D(K+1 >+0.5* (D(K+2 )-D(K+l ) ) 
?T(K)=PP(K+l)+0.5*(PP(K+2)-PP(K+l) ) 

Pi ( K ) =P l K+ 1 ) +0 . 5* t P ( K+2 ) -P ( K+l ) ) 

GO TO 205 

303 DE1(K)=DEN(K)+0.5*(DEN(K+1)-DEN<K) ) 
D1(K)*D(K)+0.5*(D(K+1)-D(K) > 

PT (K)=PP(K)+0»5*(PP( K+ 1 ) — PP ( K ) ) 
P1(K)=P(K)+0.5*(PIK+1)-P(K) ) 

GO TO 205 

218 DE1(K)=DEN(K+1 ) $ Dl ( K ) =D( K+l ) 

PT ( K ) =PP ( K + l ) 

PI ( K ) =P ( K+l ) 

GO TO 205 

219 DE1(K)=DEN(K) $ D1(K)=D(K) 

PT ( K ) =PP ( K ) 

Pi { K ) =P ( K ) 

205 CONTINUE 
DO 206 K=KB21 . KBF 1 
DE|1( K)=DEN(K) $ D1(K)=D(K) 

P T-C K ) =PP ( K ) 

206 P1;(K)=P(K> 

GO' TO 63 
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202 DE1IKB1 >=DEN(KB> 

F’T ( KB1 ) =PP (KB) 

D1 ( KB1 ) =0 ( KB > $ PI ( KB1 ) *P ( KB ) 

DO 208 K = KB , KBF 1 
IFIK.lE.KB+2) 209. 210 

209 DE 1 ( K ) =DEN ( K+2 ) $D1 ( K ) =D I K+2 ) 

'’T ( K ) =PP ( K + 2 ) 

PI ( K ) =P ( K+2 ) 

GO TO 208 

210 IF(K.lT.KB+4) 211. 212 

211 DE 1 ( K ) = DEN ( K+2 ) 

DllK)=DtK+2) 

P T ( K ) =PP ( K+2 ) 

° 1 ( K ) =P ( K+2 ) 

GO TO 208 

21 2 IFU.LE.KB+16) 213. 214 

213 DEl(K) =DEN(K+4) 

D1 ( K ) =D ( K+4 ) 

PT I K ) =PP ( K+4 ) 

P 1 ( K ) =P I K+4 ) 

GO TO 208 

214 IFIK.LE.KB+20) 215. 216 

215 L=K/2 $ IF I K.EQ .2*L)221, 222 

221 IF(K.EQ.KB;*-17) 281, 282 

281 DEl(K)=DEN(K+3)+0.5*(D£N(K+4)-DEN(K+3> ) 
Dl(K)=D(K+3)+0.5*(D(K+4)-D(K+3) ) 

PT(K)=PP(K+3)+0*5*(PP(K+4) — PP ( K + 3 ) ) 
Pl(K)=P(K+3)+0.5*(P<K+4)-P<K+3) ) 

GO TO 208 

2 82 DEI I K)=DEN I K+2 1+0.5* < DEN (K + 3 1 -DENI K+2) ) 

D1 IK)=DIK+2) +0.5*1 D (K+3 )-D (K+2 ) ) 

PT (K ) =PP( K+2 )+0.5*(PP( K+3 ) — PP { K+2 ) ) 

Pl(K)=P(K+2)+0.5*(P( K+3 ) — P I K + 2 > ) 

GO TO 208 

222 lFIK.EQ.KB+18) 284, 285 

284 DE 1 ( K ) =DEN (K+3)$Dl(K)=D(K+3) $P 1( K) =P (K+3 )$PT I K ) =PP I K+3 ) SGO TO 208 

285 DE1(K)=DEN(K+2)SD1(K)=D( K+2 )$P1 (K)=P (K+2 ) $PT I K ) =PP I K+2 ) i GO TO 208 

216 L = K/2 $ IF(K.EQ.2*L) 223, 224 

223 IFIK.EQ.KB+21) 286. 287 

286 DEl(K)=DEN(K+l)+0.5*<DEN(K+2)-DEN(K+l) ) 
l’l I K) =D< K+l )+0.5*(D(K+2 )-D< K + l ) > 

PT(K)=PP(K+l)+0.5*(PP(K+2) — PP I K+ 1 ) ) 

Pl(K)=P(K+l)+0«5*(P(K+2) — P (K+l ) ) 

GO TO 208 

287 DEI (K)=DEN(K)+0.5*( DENI K+l )-DEN(K> ) 

01(K)=D(K)+0.5*(D(K+1)-D(K) ) 

PT (K)=PP(K)+0.5*(PP(K+1) —PP I K) ) 

PI I K ) =P ( K ) +0. 5* I P I K+ 1 > — P IK)) 

GO TO 208 

224 DE 1 ( K > =DE 1 ( K + l ) 

D1 1 K ) =D I K+ 1 ) 

PT(K)=PP(K+1) 

P1(K)=P(K+1) 

208 CONTINUE 
63 CONTINUE 
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U U 


RETURN 

END 

SUBROUTINE OHO 

C THE CALCULATION OF FILM THICKNESS. DENSITY AND VISCOSITY. 

COMMON KI .KQ.KB.KBF.KF.MI .MF.KH.KKI . K2 » KKO » KKB » K02 » KKF » KF2 . 
lK3.MMI.KKH.Ki3I , KBFF , I S , N .NN , XH »U2 , PO . C2 . RC . X3 , PC » ALPH , E , P I I , 
2P1»H»HH,P.DEN,DEND,V1S,VISD,DE1,H1,PAI .VD.PSl .PSD.PS.DENX.C, 
3KW.I , J,K,XC,KB2,K31,EN,ED»PH.M.EE,SUBV.KIB.KIBF,C1.DX.S, IIS. 
4KA.S1.H0.DT ,U1 .R.Q.D.DP.DENT .C9.BET .D1 .NCT .C3.C4.U.DEO.C6.SOE 
5PT.PP.S3.S4 

DIMENSION XH( 65) ,R(40»40 ) .EE( AO) .H( 57) .HI I 57) »HH( 57) .P( 65) . 

IP 1 (57) , DENI 57) .DEND157 ) .SOE( 57 ) .PAI ( 57 > . VI S ( 57 ) . VISD ( 57 ) . 

2D (57) »PS( 57) .DX( 57) .DP (57) .VD( 57) .At 57) »S( 57) »Ql 57.65 ) . 

3DE1( 57 ) »HO(21 ) ,SUBV<65 ) , F (65 ) . C < 2 ,3 *65 ) » BET ( 40 .40 > »Dl ( 57 ) . 

4PP ( 57 ) » PT ( 57 ) 

IF(IIS.EQ.I) 18. 41 
41 I F ( I S- 2 » 10,11, 12 

10 DEN(K)=1.0+EN*P(K)/(1.0+ED*P(K) ) 

DEND(K) =EN/ ( l 1 .0+P ( K ) *ED ) **2 ) 

VIS(K)=EXP(ALPH*P(K) » 

VlSD(K)=ALPH*VIS(K) 

60 TO 25 

11 DO 13 K = KI, KO 
DEN(K)=1.0+EN*P(K)/< 1.0+ED*P(K) ) 

DEND ( K ) = EN/ ( ( 1 . 0+P ( K ) *ED ) **2 > 

VI S ( K ) = EXP ( ALPH*P ( K ) ) 

VISD(K)=ALPH*VIS(K) 

1.3_C.0NJJ_NUE_ 

GO TO 18 

12 DO 20 K = KA. KO 
P4=(P(K+1)+P(K) ) *0.5 

DEN( K )=1.0+EN*P4/ ( 1.0+ED*P4) 

QEND(K)=0.5*EN/< ( 1 .0+ED*P4 ) **2 ) 

VIS(K)=EXP( ALPH*P4) 

VISD(K) =0.5*ALPH*VIS(K) 

20 CONTINUE 
18 DO 14 K = KI ♦ KO 
DS = 0 .0 

DO 15 J=K I » KF 
15 DS=DS+Q(K, J)*P ( J) 

D(K)=-C2*DS 

H(K)=1.0+C1*0.5*(XH(K)**2)+D(K) 

IF (K.GE.KBI . AND. K.LE. KBFF ) 31, 32 

31 HH (K)=H(K)+XC*( (ABS(X3-XH(K) ) )**2-C9**2) $ GO TO 14 

32 HH ( K ) =H ( K ) 

14 CONTINUE 
25 CONTINUE 

RETURN 
END 

SUBROUTINE INTEG1 < A , B , KCT ,F ,NP , VALUE , I ERR ) 

THE CALCULATION OF INTEGRALS BY THE OVERLAPPING PARABOLA 
FORMULA. 
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COMMON K I ,KO,KB,KBF,KF,MI »MF,KH,K.KI * <2 * KK.0 * KKB * K02 * KK.F * KF2 ♦ 
1X3* MM I * K.K.H »KBI ,KBFF » I S ,N ,NN , XH ,U2 ,PO , C2 » RC » A3 » PC » ALPH , E , P I I * 
2P1 *H.HH.P.DEN»DEND.VIS,VlSD,OEl.HltPAI *VD.PSI *PSD*PS»DENX.C » 
3KW * I »J,K»XC»K.a2»KBl»EN,ED»PH,M»EE»SUBV»KIB»KIBF»Cl .DX.S* I IS* 
4KA »S1 »HO»DT. Ul, R.Q.D.DP, DENT »C9, BET. D1.NCT.C3.C4.U, DEO, C6.SOE, 
5PT »PP »S3 * S4 

DIMENSION XH(65) *R(40.40) *EE(40) *H(57) »H 1 1 5 7 ) »HH(57> *P(65) . 

1P1 ( 57) , DENI 57 ) ,DEND( 57 ) ,SOE( 57 ) *PAI ( 57 ) * VI S ( 57 ) * V I SD ( 57 ) . 

2D( 57 ) ,PS( 57 ) *DX( 57) .DPI 57) ,VDf 57) »A( 57) »S( 57) ,Q( 57,65 ) , 

3DE1I 57 ) .HO (21) , SUBV ( 65 ) ,F(65),C(2.3.65).BET(40.40)»D1(57), 

4PPI 57) ,PT( 57) 

IF (NP.LE.3) GO TO 96 

* CALCULATION OF INTERVALS OF X 

nh=np-i 

DO 10 1=1, NH 
10 DX( I ) =XH ( I + 1 ) — XH I I ) 

DO 20 1=1, NH 

* DEFINE COEFFICIENTS OF FIRST PARABOLA 
I F I I.EQ.l) GO TO 15 

C( 1.1 » I ) =-DX ( I ) **3/ I 6 * 0*DX I 1-1 )*<DX< 1-1 )+DX< I ) ) ) 

C< 1,2, I ) =DX< I)*( 3.0*DX< 1-1 )+DX( I ) ) / I 6.0*DX < I — 1 ) ) 

C I 1,3, I )=DX( I )*( 3.0*DX( I-1)+2.0*DX( I > ) / I 6 . 0* ( DX ( I-D+DXI I ) ) ) 

15 CONTINUE 

I F ( I .EQ.NH) GO TO 20 

C I 2 , 1 . 1 ) =DX ( I)*(2.0*DX< I ) + 3 • 0*DX ( I + 1 ) > / ( 6 . 0* ( DX ( I )+DX( 1 + 1) ) ) 

C( 2,2, I ) =DX ( I ) * ( DX ( I ) + 3 • 0*DX I 1 + 1 ) )/( 6.0*DX< 1+1 ) ) 

C I 2.3. I )=-DX< I )**3/ (6.0+DXI 1 + 1 ) * I DX ( I )+DX( 1 + 1) ) ) 

20 CONTINUE 

ENTRY INTEG2 

VALUE=0.0 

IF IB-A) 40,92,30 

* B IS GREATER THAN A 
30 ALIM = A 

3LIM = B 
SIGN = 1.0 
GO TO 50 

* A IS GREATER THAN B 
40 ALIM = B 

8LIM = A 
SIGN =-1.0 
50 NH=NP-1 

IF(KCT.EQ.l) 125, 123 

* CALCULATION OF INTEGRAL OVER SUBINTERVAL 
123 DO 80 1=1, NH 

SUBV I I ) =0.0 

IF IXHI I ) .EQ.ALIM) SUBVI I ) =C 12 , 1 , 1 ) *F ID +C I 2 , 2 , 1 ) *F ( I + 1 ) 

X+C (2,3,1 )*F! 1+2) 

IF IXHI 1+1 ) .EQ.BLIM) SUBVI I )=C( 1 , 1 , I ) *F I I -1 ) +C I 1 ,2 , I ) *F I I ) 

X+CI 1.3, I >*F< 1+1) 

IF IXHI I ). GT. ALIM. AND. XH< I + D.LT.BLIM) SUBVI I ) =0. 5* I C 1 1 , 1 . 1 ) 

X*F I I-1)+(C(1,2,I)+C(2,1,I))*F(I)+(CI1,3,I)+C(2,2,I))*F(I+1)+ 
XC(2,3*I )*F( 1+2) ) 

80 VALUE=VALUE+SUBV( I ) 

VALUE=SIGN*VALUE 
GO TO 92 
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r> n 


125 DO 110 1=1. NH 
SUBV( I ) =0.0 

IF(XH( I J.EQ.ALIM) 111. 110 
111 IF(I.EQ.NH) 113. 114 

113 SUBV ( I)=C(1»1»I)*F( I-l)+C (1.2.1 )*F( I )+C (1.3.1 ) *F ( I +1 ) 

GO TO 120 

114 I F ( I «GE « 2 ) 115, 116 

115 SUBV ( I )=0.5*(C( 1.1, I )*F( 1-1 J + (C( 1,2, I )+C(2.1.I ) )*F( I ) + (C( 1 .3. 1 ) 
X+C (2,2.1) )*FII+1)+C(2»3»I ) *F ( 1+2) ) 

GO TO 120 

116 SUBV ( I)=C(2»l»I)*F(I)+C(2»2»I)*F(I+l)+C(2.3»I)*F(I+2) 

120 VAU/E = SIGN*SUBV( I) $ GO TO 92 

110 CONTINUE 

* SET ERROR PARAMETER FOR TOO FEW POINTS 
92 IERR = 0 

RETURN 

* SET ERROR PARAMETER FOR NORMAL RETURN 

96 IERR = 1 
RETURN 

* SET ERROR PARAMETER FOR A AND/OR B OUT OF RANGE OF TABLE 

97 IERR = 2 
RETURN 
END 

SUBROUTINE INTEG 

THE CALCULATION OF (IX), THE INTEGRAL OF 
SQUEEZING TERM. 

COMMON XI .XO.XB.XBF.XF.MI «MF .XH.XXI , K2 . XXO » XXB , X02 .XXF.XF2, 

1X3, MM I .XXH.XBI , XBFF , I S , N , NN ,XH , U2 , PO , C2 , RC » X3 , PC , ALPH , E , P I I , 

2P-1-. H., HH-.P , D E N , D E N D , V I S . V IS D.DEl.Hl.PAI.VD.PSI.PSD.PS.DENX.C. 

3XW.I ,J,X,XC»XB2,XB1»EN,ED»PH, M7E EYS U B VTXTBtK-I- B F ,-C-l rD-X.-S-.-I -I-S-. 

4XA»S1»H0,DT»'J1»R»Q»D.DP»DENT»C9.BET»D1.NCT»C3.C4*U»DE0,C6*S0E» 
5PT.PP.S3.S4 

DIMENSION XH(65) ,R( 40,40 ) ,EE<40 ) .H( 57) »H1 ( 57) »HH( 57 ) »P( 65 ) » 
1PK57) .DENI 57) ,DEND(57) ,SOE<57> ,PAI(57)»VIS(57),VISD(57)» 

2D (57) .PS (57 ) ,DX< 57) ,DP( 57) .VD( 57 ) .At 57) »S< 57) ,Q< 57.65 ) , 

3DEK57.) »HO (21) .SUBV (65) ,F(65)»C(2,3»65>»BET(40,40),D1(57)» 

4PP ( 57 ) ,PT ( 57 ) 

IF(M.EQ.I) 81, 82 

81 DO 83 X=XI » XO 
83 PA I ( X ) =0 • 0 $ GO TO 86 

82 DO 60 X=X I , XO 

SD=(HO(M)*DEN(X)-HO(M-l )*DEl (X) )/DT 
202 IF(X.GE.XBI. AND. X.LT. XBFF) 61, 62 

61 DENT= ( DEN ( X ) -DE 1 ( X ) ) /DT 
DENX= ( DEN < X ) -DEN ( X-l ) ) /DX ( X- 1 ) *0 . 5* ( S4+S3 ) 

SE=XC*HO(M>*( ( ABS ( X3-XH ! X ) > ) **2-C9**2 ) * ( DENT+DENX ) 

SOE ( X ) =SD+SE 
GO TO 60 

62 SOE ( K ) =SD 
60 CONTINUE 

PA I ( XO) =0.0 
X01=X0-1 
DO 63 X=XI » XOl 
J=XO-X 
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63 PA I < J ) = PAI ( J+l )+0.5*DX< J ) * < SOE < J ) +SOE { J+l) ) 

I F ( I IS.EQ.l) GO TO 86 
DO 101 K=KA , KO 

00 101 J=KA » KO 

SSS=0.0 $ NN*K-KA+1 $ II=J-KA+1 
8ET(NN, I I )=0.0 
IF(K.EQ.KO) GO TO 161 
IF(K.GE.KBI .AND.K.LT.KBFF) 162,155 

162 IF(J.GE.K) 154,155 

154 IF(K.EQ.KBFF-1 ) 163,164 

163 I F ( J.EQ.K) 165, 155 

165 BET ( NN , 1 1 )=XC*HO<M)*< ( ABS ( X3“XH ( J ) ) ) **2~C9**2 ) *DEND ( J ) * 

X ( 0 . 5* ( S4+S3 ) *DX ( J ) /DX( J-l ) +DX ( J ) /DT ) *0. 5 
GO TO 155 

164 BET(NN,II )=XC*HO(M)*( ( ABS(X3-XH(J) ) )**2-C9**2) 

X*DENO( J)*(0.5*< S4+S3 )/DX( J-l J+1.0/DT)*(XH( J+l )-XH( J-l ) ) * 0 • 5 

155 IF(J-K) 101, 107, 108 

107 BET ( NN, I I )*0.5*DX( J )*DEND( J ) *H0 { M ) /DT+BET ( NN , I I ) 

GO TO 101 

108 BET ( NN, 1 I )=0.5*DEND( J ) *H0 < M ) * { XH ( J + l ) -XH C J-l ) ) /DT + BET (NN, I I ) 

GO TO 101 

161 BET(KO-KA+1,II)=0.0 

101 CONTINUE 
86 CONTINUE 
RETURN 

203 FORMAT ( 4H FF=,E12.5,4H F2=»E12.5,4H F3=.E12.5»4H F4=,E12.5, 

X4H F5=»E12#5) 

204 FORMAT ( 6HH0 (M)=»E12.5,4H C3= » E 12 . 5 »6HXH ( K ) = » E 1 2 . 5 . 7HDEN ( K ) =E 1 2 . 5 ) 

205 FORMAT ( 8HH0 ( M— 1)=,E12»5, 4H DD*,E12.5,6H DDD= , E 12 . 5 , 7HDE 1 ( K ) = ♦ 
XE12.5) 

210 FORMAT (50X»*KERNEl = *//(5(2I4»E14,5) ) ) 

211 FORMAT <50X»*bEND(K)=*//(6( I4»1X,E15.7))! 

212 FORMAT < 50X»*DX(K)= + // <6( I4,1X,E15.7))) 

END 

FUNCTION PHI < SR ) 

C THE CALCULATION OF THE COEFFICIENTS OF EQUAT I ON (. 54 ) OF PART II. 

COMMON KI ,K0,KB,KBF,KF,MI ,MF,KH,KKI , K2 » KKO , KKB , K02 , KKF , KF2 , 

1 K3 » MM I , KKH , KB I , KBFF , I S , N » NN , XH , U2 , PO * C2 » RC , X3 , PC , ALPH , E , P I I , 
2P1»H,HH,P,DEN»DEND,VIS,VISD,DE1»H1 *PAI , VO, PS I , PSD , PS , DENX ,C , 

3KW , I »J,K,XC»KB2.K81»EN,ED,PH,M,EE,SUBV,KIB,KIBF,C1 ,DX,S, I IS, 

4KA *S1 »H0, DT , Ul, R, Q, D,DP, DENT, C9*8ET,D1,NCT,C3»C4,U, DEO, C6» SOE, 

5PT ,PP ,S3,S4 

DIMENSION XH(65) ,R(40,40) » EE ( 40 ) * H ( 57) » H 1 ( 5 7 ) * HH (57) , P l 65 ) » 

1P1 ( 57) , DEN ( 57! ,DENO( 57) ,SOE( 57 ),PAI(57)»VIS(57),VISD(57), 
2D(57),PS(57),DX(57),DP(57),VD(57)»A(57)*S(57),Q(57,65), 

3DE1 ( 57 ) »H0( 21 ) ,SUBV( 65 ),F(65),C(2,3»65)»BET(40,40)»D1(57), 

4PP( 57) ,PT( 57 ) 

0K0=( 1.0-1. O/EXP (ALPH) )/ALPH 
Cl=C3/SR $ C2=C1/PI I % XC=C6/SR 

1 I S= 1 $ CALL DHD $ CALL INTEG 
IF(M.EQ.I) 3, 4 

3 DO 5 K= K I * KO 

IF(K.GE.KBI. AND. K.LE. KBFF) 20, 21 
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20 STO= I S4-S3 ) *XC*< X3*XH< XBFF ) -0 . 5*XH ( XBFF ) **2-X3*XH( K ) 
X+0.5*XHIX)**2) $ GO TO 5 

21 STO=0 • 0 

5 VD(K)=( (S4+S3)*0«5*(H(X)-DE0/DENIX) > -STO ) / I HH ( X ) ** 3 ) 

CALL INTEG2(XH(KI ) ♦ XH < <0 ) ♦ 2 » VD, XO » VALUE , I ER ) 

S1=VALUE 

PHI=48.0*U/SR**2-(QX0/S1 ) 

GO TO 10 

A DO 2 X<=X I » XO 

SIX)=PAIIX)/!HH(X)**3*DEN(X) ) 

IFIX.GE.XBI. AND. X.LE. XBFF) 22 . 23 

22 ST0=IS4-S3>*XC*( X3*XH(KBFF)-0.5*XH(XBFF)**2-X3*XH< X) 
X+0.5*XH(X)**2) S GO TO 2 

23 ST0=0.0 

2 VO(X)=( ( S4+S3 ) *0 • 5* ( H ( X ) -DEO /DEN ( X ) )-ST0)/(HH(X)**3) 

CALL INTEG2(XH(KI ) »XH(XO) »2, VO, XO »VALUE , I ER ) 

Sl=VALUE 

CALL INTEG2 ( XH ( X I ) ,XH<XO) » 2 , S , XO , V ALUE » I ER ) 

S2=VALUE 

PHI=QXO*SR**3-A8.0*U*(S1*SR-S2 ) 

10 RETURN 
END 

SUBROUTINE NRONE < X ,CONV ,NI T , DELXO , XO»XMAX»XMIN) 

COMMON XI ,X0»XB»XBF»XF,MI ,MF,XH,XXI » X2 , XXO » XXB , X02 » XXF , XF 2 , 

1X3 » MM I » XXH » XB I , XBFF , I S ,N »NN , XH , U2 , PO , C2 , RC , X3 , PC , ALPH , E , P I I , 
2P1*H»HH»P»DEN»DEND,VIS,VISD,DE1,H1 ,PAI ,VD,PSI ,PSD , PS ,DENX ,C , 
3XW , I ,J,X,XC,XB2»XB1,EN,ED,PH,M,EE,SUBV,XIB,XIBF,C1,DX,S, IIS, 
4XA,S1,H0,DT,U1 ,R , Q , D ,Dp , DENT , C 9 »BET , D1 ,NCT , C3 ,CA , U ,DEO , C6 , SOE , 

5 P T,P PtS3 

DIMENSION XH(65) ,R(AO»AO) ,EE(AO) , H C 5 7) »H1( 57) ,HH( 57) , P ( 65 ) * 
1P1( 57) .DENI 57) ,DEND< 57) , SOE I 57) ,PAI (57) , VI SI 57) ,VISD(57) , 
2DI57) ,PS(57) ,DX(57) ,DP(57) ,VD(57) ,A(57) ,S(57) ,Q(57,65! , 

3DE1I 57) , HOI 21) ,SUBV ( 65 ) , F I 65 ) ,C I 2 , 3 ,65 ) , BET I AO , 40 ) ,D1 1 57 ) , 

APPI 57) ,PTI 57 ) 

X=XO $ I T= 1 $ NR=5$ NW=6$ DELX=DELXO 
WRITEI6.210) C1,C2,C9,XC,X3 
8 r C=PH I < X ) 

IFIM.EGUl) GO TO AO 
I F ( IT-1) 10, 10, 15 

10 Xl=X $ X=X+DELX $ GO TO 11 

15 FC2=FC $ DFC=(FC2-FC1)/(X2-X1) 

DELX=-FC/DFC 

IF(ABSIFC)-CONV) 25, 25, 20 
20 Xl=X $ X=X+DELX 

11 FC1=FC S X2=X $ IF(IT-NIT) 22, 22, 25 
22 I T=I T+l $ WRITE! NW,7 ) IT, FC, X, DELX 

IF(X.GE.XMAX) GO TO 25 
IFIX.LT.XMIN) GO TO 25 
GO TO 8 
AO FC=PHI (X) 

TDC=-96.0*U/X**3 
DEX=~ FC/TDC 

IFIABS(FC)-CONV) 25, 25, A2 
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42 IF(IT.LE.3> 50*51 

50 FA= I T $ GO TO 52 

51 FA=3.0 

52 0ELX=FA*D£X/3.0 

X=X+DELX $ IF(IT-NIT) 43, 43 » 25 

43 I T = I T+l $ WR I TE ( NW ,7 ) IT»FC,X,DELX 
IF(X.GT.XMAX) GO TO 25 
IF(X.LT.XMIN) GO TO 25 

GO TO 40 

25 WRITE(NW,7) IT, FC , X, DELX 
RETURN 

210 FORMAT ( 5E12 • 5 ) 

7 FORMAT ( 5H IT=, 15, 4H FC=» E12.5, 5H X= » E12.5, 
X5H SX=, E12.5) 

END 

' END OF RECORD 
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